**Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI**

Joao L. A. Carvalho1 and Krishna S. Nayak2

<sup>1</sup>*Universidade de Brasília* <sup>2</sup>*University of Southern California* <sup>1</sup>*Brazil* <sup>2</sup>*USA*

## **1. Introduction**

Aortic stenosis consists in a narrowing or incomplete opening of the aortic valve. This typically alters the blood flow, causing turbulent and/or complex flow jets. Such jets display peak velocities considerably higher than those of normal flow, and a much broader range of flow velocities. Another form of aortic valve disease is aortic insufficiency, or regurgitation. This condition occurs when the aortic valve fails to close completely. This is also known as "leaky valve", as flow jets in the reverse direction are observed when no flow should occur. The visualization and quantitation of cardiovascular blood flow is an important component of the assessment of aortic valve disease. For example, peak velocity measurements in flow jets are used to estimate pressure drop, which is an indicator of the hemodynamic load of a stenosis (Tsai et al., 1999).

#### **1.1 Doppler ultrasound**

The current non-invasive gold standard for flow quantitation is Doppler ultrasound. The ultrasound equipment is relatively inexpensive, small, and portable. Measurements are typically obtained in real-time, with excellent temporal resolution. The most popular techniques for ultrasound flow assessment are color Doppler and spectral Doppler.

Evaluation by ultrasound is impossible when there is air, bone, or surgical scar in the ultrasound path. Examination by ultrasound in obese patients is difficult, as the overlying adipose tissue (fat) scatters the sound waves. Doppler flow measurements may be inaccurate when the ultrasound beam cannot be properly aligned with the vessel axis, requiring measured velocities to be "angle-corrected" by the operator. Peak-velocity overestimation on the order of 18–40% have been reported in the literature (Hoskins, 1996; Winkler & Wu, 1995), usually due to spectral broadening at large insonation angles, and to Doppler gain settings.

#### **1.2 Magnetic resonance imaging**

Magnetic resonance imaging (MRI) is potentially the most appropriate technique for addressing all aspects of cardiovascular disease examination, which includes assessing myocardial function and perfusion, as well as visualizing and measuring blood flow. 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

region of interest (e.g., mitral valve) moves during the cardiac cycle. Other problems include

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 5

FVE has also been accelerated by simply neglecting spatial encoding along one of the spatial dimensions (Feinberg et al., 1985; Hennig et al., 1988), or by acquiring velocity images with no spatial encoding other than slice selection (Galea et al., 2002). In these techniques, the velocity measurement is a projection of all signal along a line or a plane in 3D space, respectively. As a consequence, both methods have dynamic range issues, as the signal of flowing blood has to be distinguished from all the background signal from static tissue observed along the projection. Furthermore, these approaches are unable to resolve different sources of flow that

This chapter introduces spiral FVE, a novel method for MR flow quantitation that addresses the limitations discussed above. The proposed method provides fully-localized time-velocity distribution measurements, in a single acquisition, that is one short breath-hold long (approximately 10 seconds). Spiral FVE uses conventional slice-selective excitation (Bernstein et al., 2004; Nishimura, 2010), which excites (selects) a thin slice of the body to be imaged. The two-dimensional plane defined by this slice is imaged using spiral acquisitions (Ahn et al., 1986; Meyer et al., 1992), which encode both spatial dimensions simultaneously. Therefore, no spatial encoding is neglected, and measurements are fully localized in 3D space. 2D-resolved spatial encoding allows for easy localization of the region of interest, and the ability to resolve multiple sources of through-plane flow in the imaged field-of-view, without requiring static tissue suppression. Scan-plane prescription is performed using classic protocols, which is

considerably less laborious than the beam-placing process used in real-time FVE.

Without acceleration, spiral FVE presents some limitations: (1) insufficient velocity field-of-view (the maximum range of velocities allowed without aliasing); (2) low in-plane spatial resolution, which limits the ability of spatially localizing the flow; (3) long readouts, which causes spatial blurring at 3 T, due to off-resonance effects (Noll, Meyer, Pauly, Nishimura & Macovski, 1991); and (4) moderate temporal resolution, which may blur certain features of the flow waveform. We address these limitations using the following acceleration techniques: variable-density spirals (Tsai & Nishimura, 2000), partial Fourier reconstruction (Noll, Nishimura & Macovski, 1991), and temporal acceleration (Madore et al., 1999; Tsao, 2002). By combining these techniques, we achieve a total 18-fold acceleration in

MRI is a modality uniquely capable of imaging all aspects of heart 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

The main component of a MRI scanner is a strong magnetic field, called the *B*<sup>0</sup> 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,

the large voxel size and low temporal resolution.

may co-exist in the projected line or plane.

**1.3 Chapter outline**

spiral FVE.

**2. MR flow imaging**

**2.1 Basic principles of MRI**

frequency range, which are harmless to the patient.

system. Magnetic resonance (MR) measurements are also less operator-dependent than those of Doppler ultrasound, and the true direction of flow can generally be precisely measured. The main MR techniques for measuring flow are phase contrast and Fourier velocity encoding. These techniques are introduced below, and will be discussed in further detail in section 2.

#### **1.2.1 Phase contrast**

Phase contrast (O'Donnell, 1985) is a technique in which a bipolar gradient (see section 2) aligned with the axis of flow is used to obtain a velocity measurement associated with each pixel (or "voxel") of the image. In practice, two acquisitions are used, and the first moment of the bipolar gradient is varied between measurements. The velocity estimate is obtained from the phase difference between the images obtained in each acquisition.

Phase contrast can be combined with dynamic (cine) MRI (Glover & Pelc, 1988), in which short acquisitions and some form of cardiac synchronization are used to produce images throughout the cardiac cycle. The combined technique (cine phase contrast) can depict motion and flow throughout the cardiac cycle (Nayler et al., 1986). Alternatively, phase contrast can be combined with real-time MRI (Holsinger et al., 1990; Riederer et al., 1988), in which the encodings are applied sequentially, periodically, and continuously. In real-time MRI, images are formed by sliding a window along the acquired data and reconstructing an image for each position of the window. The display of real-time phase contrast data is typically implemented as color overlay of flow information (phase difference) over the anatomical (magnitude) image, which is called real-time color flow (Nayak et al., 2000; Riederer et al., 1991).

In phase contrast, data inconsistency, partial volume effects, and intravoxel phase dispersion can lead to peak velocity underestimation (Clarke et al., 1996; Tang et al., 1993). Partial voluming is particularly problematic when flow is highly localized and/or turbulent. When a large voxel size is adopted to measure the flow rate, not only may moving spins and stationary spins coexist in a voxel, but also the velocity distribution of spins within a voxel may spread over a wide range of velocities. This results in signal loss, distortion and erroneous velocity estimates. As a result, phase contrast imaging can not provide accurate peak velocity measurements in turbulent and/or complex flow jets. Such jets are commonly observed in narrowed vessels, and in valves presenting stenotis and/or regurgitation.

#### **1.2.2 Fourier velocity encoding**

The limitations mentioned above can be overcome using Fourier velocity encoding (FVE) (Moran, 1982). FVE can be considered the MR equivalent to spectral Doppler. In this technique, the full spectrum of velocities within each voxel is measured by phase-encoding the velocity information in Fourier domain. Therefore, FVE is robust to partial voluming, and flow measurements from low spatial resolution images are still accurate (Tsai et al., 1999). FVE shows satisfactory agreement with Doppler ultrasound (Mohiaddin et al., 1997). However, it is typically not used clinically, because the acquisition time required by this technique is in principle considerably longer than that of phase contrast.

Different approaches to accelerating FVE have been proposed. One example is the use of two-dimensional cylindrical excitation to restrict the field-of-view to a beam that can be imaged without phase encoding (Dumoulin et al., 1991). This approach makes it possible to perform spatial encoding and velocity encoding simultaneously, and in a single pulse repetition time (TR) (DiCarlo et al., 2005; Hu et al., 1993; Irarrazabal et al., 1993; Macgowan et al., 2005). This allows FVE measurements to be obtained in real-time. However, real-time FVE has problems related to the precise placement of the imaging beam, especially when the region of interest (e.g., mitral valve) moves during the cardiac cycle. Other problems include the large voxel size and low temporal resolution.

FVE has also been accelerated by simply neglecting spatial encoding along one of the spatial dimensions (Feinberg et al., 1985; Hennig et al., 1988), or by acquiring velocity images with no spatial encoding other than slice selection (Galea et al., 2002). In these techniques, the velocity measurement is a projection of all signal along a line or a plane in 3D space, respectively. As a consequence, both methods have dynamic range issues, as the signal of flowing blood has to be distinguished from all the background signal from static tissue observed along the projection. Furthermore, these approaches are unable to resolve different sources of flow that may co-exist in the projected line or plane.

#### **1.3 Chapter outline**

2 Will-be-set-by-IN-TECH

system. Magnetic resonance (MR) measurements are also less operator-dependent than those of Doppler ultrasound, and the true direction of flow can generally be precisely measured. The main MR techniques for measuring flow are phase contrast and Fourier velocity encoding. These techniques are introduced below, and will be discussed in further detail in section 2.

Phase contrast (O'Donnell, 1985) is a technique in which a bipolar gradient (see section 2) aligned with the axis of flow is used to obtain a velocity measurement associated with each pixel (or "voxel") of the image. In practice, two acquisitions are used, and the first moment of the bipolar gradient is varied between measurements. The velocity estimate is obtained from

Phase contrast can be combined with dynamic (cine) MRI (Glover & Pelc, 1988), in which short acquisitions and some form of cardiac synchronization are used to produce images throughout the cardiac cycle. The combined technique (cine phase contrast) can depict motion and flow throughout the cardiac cycle (Nayler et al., 1986). Alternatively, phase contrast can be combined with real-time MRI (Holsinger et al., 1990; Riederer et al., 1988), in which the encodings are applied sequentially, periodically, and continuously. In real-time MRI, images are formed by sliding a window along the acquired data and reconstructing an image for each position of the window. The display of real-time phase contrast data is typically implemented as color overlay of flow information (phase difference) over the anatomical (magnitude)

image, which is called real-time color flow (Nayak et al., 2000; Riederer et al., 1991).

narrowed vessels, and in valves presenting stenotis and/or regurgitation.

principle considerably longer than that of phase contrast.

In phase contrast, data inconsistency, partial volume effects, and intravoxel phase dispersion can lead to peak velocity underestimation (Clarke et al., 1996; Tang et al., 1993). Partial voluming is particularly problematic when flow is highly localized and/or turbulent. When a large voxel size is adopted to measure the flow rate, not only may moving spins and stationary spins coexist in a voxel, but also the velocity distribution of spins within a voxel may spread over a wide range of velocities. This results in signal loss, distortion and erroneous velocity estimates. As a result, phase contrast imaging can not provide accurate peak velocity measurements in turbulent and/or complex flow jets. Such jets are commonly observed in

The limitations mentioned above can be overcome using Fourier velocity encoding (FVE) (Moran, 1982). FVE can be considered the MR equivalent to spectral Doppler. In this technique, the full spectrum of velocities within each voxel is measured by phase-encoding the velocity information in Fourier domain. Therefore, FVE is robust to partial voluming, and flow measurements from low spatial resolution images are still accurate (Tsai et al., 1999). FVE shows satisfactory agreement with Doppler ultrasound (Mohiaddin et al., 1997). However, it is typically not used clinically, because the acquisition time required by this technique is in

Different approaches to accelerating FVE have been proposed. One example is the use of two-dimensional cylindrical excitation to restrict the field-of-view to a beam that can be imaged without phase encoding (Dumoulin et al., 1991). This approach makes it possible to perform spatial encoding and velocity encoding simultaneously, and in a single pulse repetition time (TR) (DiCarlo et al., 2005; Hu et al., 1993; Irarrazabal et al., 1993; Macgowan et al., 2005). This allows FVE measurements to be obtained in real-time. However, real-time FVE has problems related to the precise placement of the imaging beam, especially when the

the phase difference between the images obtained in each acquisition.

**1.2.1 Phase contrast**

**1.2.2 Fourier velocity encoding**

This chapter introduces spiral FVE, a novel method for MR flow quantitation that addresses the limitations discussed above. The proposed method provides fully-localized time-velocity distribution measurements, in a single acquisition, that is one short breath-hold long (approximately 10 seconds). Spiral FVE uses conventional slice-selective excitation (Bernstein et al., 2004; Nishimura, 2010), which excites (selects) a thin slice of the body to be imaged. The two-dimensional plane defined by this slice is imaged using spiral acquisitions (Ahn et al., 1986; Meyer et al., 1992), which encode both spatial dimensions simultaneously. Therefore, no spatial encoding is neglected, and measurements are fully localized in 3D space. 2D-resolved spatial encoding allows for easy localization of the region of interest, and the ability to resolve multiple sources of through-plane flow in the imaged field-of-view, without requiring static tissue suppression. Scan-plane prescription is performed using classic protocols, which is considerably less laborious than the beam-placing process used in real-time FVE.

Without acceleration, spiral FVE presents some limitations: (1) insufficient velocity field-of-view (the maximum range of velocities allowed without aliasing); (2) low in-plane spatial resolution, which limits the ability of spatially localizing the flow; (3) long readouts, which causes spatial blurring at 3 T, due to off-resonance effects (Noll, Meyer, Pauly, Nishimura & Macovski, 1991); and (4) moderate temporal resolution, which may blur certain features of the flow waveform. We address these limitations using the following acceleration techniques: variable-density spirals (Tsai & Nishimura, 2000), partial Fourier reconstruction (Noll, Nishimura & Macovski, 1991), and temporal acceleration (Madore et al., 1999; Tsao, 2002). By combining these techniques, we achieve a total 18-fold acceleration in spiral FVE.

#### **2. MR flow imaging**

#### **2.1 Basic principles of MRI**

MRI is a modality uniquely capable of imaging all aspects of heart 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*<sup>0</sup> 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,

When the readout gradients are played, the acquired signal at a particular time instant corresponds to the sum of different sinusoidal signals generated by spins located at different regions of the body, each rotating at different frequencies corresponding to their spatial locations. If an axial slice is being acquired, for example, the demodulated signal value is equivalent to a sample of the Fourier transform *M*(*kx*, *ky*) of the cross-sectional image *m*(*x*, *y*). In this case, by changing the amplitudes of *Gx* and *Gy* during acquisition, one may acquire different samples of *M*(*kx*, *ky*). In fact, by playing *Gx* and/or *Gy*, one can move along the *kx*-*ky* plane (which is known in MRI as k-space), collecting samples of *M*(*kx*, *ky*). When enough samples of *M*(*kx*, *ky*) 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 *kx*-*ky* 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

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 7

enough to cover the entire object, overlap in spatial domain is observed (aliasing).

called pulse repetition time, or TR. **a**

**b**

and (b) spiral acquisitions.

Gz

RF

Gz

RF

Gx

Gy

Gx

Gy

Because signal amplitude is lost as the net magnetization realigns with *B*<sup>0</sup> (this is called relaxation), multiple acquisitions (excitation + readout) may be needed in order to cover the entire k-space. Different 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 *kx*-*ky* 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 fashion in which RF pulses and gradients are played is called pulse sequence. The time between acquisitions is

Fig. 1. Timing diagram (left) and corresponding k-space trajectories (right) for (a) 2DFT,

kx

kx

ky

ky

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*<sup>0</sup> 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 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*<sup>0</sup> field, so MRI scanners with stronger magnetic fields (e.g., 3 Tesla) provide higher signal-to-noise ratio (SNR).

Another important component of the scanner are the gradient coils. There are typically three gradient coils (*Gx*, *Gy*, and *Gz*), that produce an intentional perturbation in the *B*<sup>0</sup> field when turned on ("played"). This perturbation varies linearly along each spatial direction (*x*, *y* and *z*), such that no perturbation is perceived 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*<sup>0</sup> 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, *Gx*, *Gy*, and *Gz* 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 a 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*<sup>0</sup> 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 (*Mxy*).

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 an 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, *Mxy* 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) , *Gz* is played during excitation (for slice-selection), and *Gx* and *Gy* 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.

4 Will-be-set-by-IN-TECH

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*<sup>0</sup> 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 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*<sup>0</sup> field, so MRI scanners with stronger magnetic fields (e.g., 3 Tesla) provide higher signal-to-noise ratio (SNR). Another important component of the scanner are the gradient coils. There are typically three gradient coils (*Gx*, *Gy*, and *Gz*), that produce an intentional perturbation in the *B*<sup>0</sup> field when turned on ("played"). This perturbation varies linearly along each spatial direction (*x*, *y* and *z*), such that no perturbation is perceived 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*<sup>0</sup> 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, *Gx*, *Gy*, and *Gz* are used to frequency-encode (or phase-encode) spatial position along the *x*, *y* and *z* directions,

The final major component of the MR scanner is the radio-frequency (RF) coil. This is used to transmit a 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*<sup>0</sup> 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

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 an 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

When the RF pulse is turned off, *Mxy* 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) , *Gz* is played during excitation (for slice-selection), and *Gx* and *Gy* are played during acquisition. These can be switched, for acquiring sagittal or coronal images, or all three gradients may be used during both excitation

respectively.

combination.

will now have a component in the *x*-*y* plane (*Mxy*).

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 different sinusoidal signals generated by spins located at different regions of the body, each rotating at different frequencies corresponding to their spatial locations. If an axial slice is being acquired, for example, the demodulated signal value is equivalent to a sample of the Fourier transform *M*(*kx*, *ky*) of the cross-sectional image *m*(*x*, *y*). In this case, by changing the amplitudes of *Gx* and *Gy* during acquisition, one may acquire different samples of *M*(*kx*, *ky*). In fact, by playing *Gx* and/or *Gy*, one can move along the *kx*-*ky* plane (which is known in MRI as k-space), collecting samples of *M*(*kx*, *ky*). When enough samples of *M*(*kx*, *ky*) 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 *kx*-*ky* 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 is observed (aliasing).

Because signal amplitude is lost as the net magnetization realigns with *B*<sup>0</sup> (this is called relaxation), multiple acquisitions (excitation + readout) may be needed in order to cover the entire k-space. Different 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 *kx*-*ky* 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 fashion in which RF pulses and gradients are played is called pulse sequence. The time between acquisitions is called pulse repetition time, or TR.

Fig. 1. Timing diagram (left) and corresponding k-space trajectories (right) for (a) 2DFT, and (b) spiral acquisitions.

are moving.

**2.3.1 Phase contrast**

**2.3.2 Fourier velocity encoding**

moment of *G <sup>r</sup>*(*t*):

velocity times the first moment of the gradient waveform along the direction in which they

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 9

For spins moving along the *r*-axis with a constant velocity *v*, and initial position*r*0, we can

*<sup>G</sup> <sup>r</sup>*(*t*) *dt* <sup>+</sup> *<sup>γ</sup><sup>v</sup>* ·

where *M* <sup>0</sup> and *M* <sup>1</sup> are the zeroth and first moments of the*r*-gradient waveform at echo time, respectively. Thus, if a gradient with null zeroth moment is used (e.g., a bipolar gradient,

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

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:

> *<sup>v</sup>*(*x*, *<sup>y</sup>*) = *<sup>φ</sup>a*(*x*, *<sup>y</sup>*) <sup>−</sup> *<sup>φ</sup>b*(*x*, *<sup>y</sup>*) *γ*(*M<sup>a</sup>*

While phase contrast provides a single velocity measurement associated with each voxel, Fourier velocity encoding (FVE) (Moran, 1982) provides a velocity histogram for each spatial

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*kv* is the velocity frequency variable associated with *v*, and is proportional to the first

*kv* <sup>=</sup> *<sup>γ</sup>* 2*π*

where *φa*(*x*, *y*) and *φb*(*x*, *y*) are the phase images acquired in each acquisition, and *M<sup>a</sup>*

location, which is a measurement of the velocity distribution within each voxel.

*φ*(*r*,*v*, *t*) = 2*π*(

<sup>1</sup> <sup>−</sup> *<sup>M</sup><sup>b</sup>*

*<sup>G</sup> <sup>r</sup>*(*t*) · (*<sup>r</sup>*<sup>0</sup> <sup>+</sup>*vt*) *dt* (9)

*G <sup>r</sup>*(*t*)*t dt* (10)

<sup>1</sup>) , (12)

*kr* ·*<sup>r</sup>* <sup>+</sup>*kv* ·*<sup>v</sup>*), (13)

*M* 1. (14)

<sup>1</sup> and *<sup>M</sup><sup>b</sup>* 1

 *t*<sup>0</sup> 0

<sup>=</sup> *<sup>γ</sup><sup>r</sup>*<sup>0</sup> · *<sup>M</sup>* <sup>0</sup> <sup>+</sup> *<sup>γ</sup><sup>v</sup>* · *<sup>M</sup>* 1, (11)

write*r*(*t*) = *r*<sup>0</sup> +*vt*. Rewriting equation 6, for *t* = *t*0:

*φ* = *γ*

= *γr*<sup>0</sup> ·

 *t*<sup>0</sup> 0

> *t*<sup>0</sup> 0

aligned with *<sup>v</sup>*), the phase accrued for a constant velocity spin is *<sup>φ</sup>* <sup>=</sup> *<sup>γ</sup><sup>v</sup>* · *<sup>M</sup>* 1.

erroneous interpretation of the local velocity (Rebergen et al., 1993).

are the first moment of the bipolar gradients used in each acquisition.

#### **2.2 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*(*kx*, *ky*) of the cross-sectional image *m*(*x*, *y*):

$$M(k\_{\lambda}, k\_{y}) = \int\_{\mathcal{X}} \int\_{\mathcal{Y}} m(\mathbf{x}, y) \, e^{-j2\pi(k\_{\lambda}\mathbf{x} + k\_{y}y)} \, d\mathbf{x} \, dy. \tag{1}$$

The Fourier coordinates *kx* and *ky* vary with time, according to the zeroth moment of the readout gradients *Gx* and *Gy*:

$$k\_{\mathbf{x}}(t) = \frac{\gamma}{2\pi} \int\_{0}^{t} G\_{\mathbf{x}}(\tau) \, d\tau \tag{2}$$

$$k\_{\mathcal{Y}}(t) = \frac{\gamma}{2\pi} \int\_0^t \mathbb{G}\_{\mathcal{Y}}(\tau) \, d\tau. \tag{3}$$

These equations explain how the gradients can be used to "move" along k-space, as discussed on section 2.1. This formalism can be generalized for any combination of the gradients *Gx*, *Gy* and *Gz* as:

$$M(\vec{k}\_r) = \int\_{\vec{r}} m(\vec{r}) \cdot e^{-j2\pi \vec{k}\_r \cdot \vec{r}} \, d\vec{r} \tag{4}$$

$$
\vec{k}\_r(t) = \frac{\gamma}{2\pi} \int\_0^t \vec{G}\_r(\tau) \, d\tau,\tag{5}
$$

where *G <sup>r</sup>* is the oblique gradient resulting from the combination of the *Gx*, *Gy* and *Gz* gradients, and*r* is its corresponding axis, along which the linear variation in magnetic field intensity is perceived.

Given a spatial position function*r*(*t*) and a magnetic field gradient *G <sup>r</sup>*(*t*), the magnetization phase is:

$$\Phi(\vec{r},t) = \gamma \int\_0^t \vec{G}\_r(\tau) \cdot \vec{r}(\tau) \,d\tau,\tag{6}$$

For static spins,*r*(*t*) is constant (*r*), and this becomes:

$$
\phi = \gamma \vec{r} \cdot \int\_0^t \vec{G}\_r(\tau) \, d\tau \tag{7}
$$

$$
\vec{r} = 2\pi \vec{k}\_r \cdot \vec{r}\_r \tag{8}
$$

as in the exponential in equation 4.

#### **2.3 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 1950's. However, clinical applications of MR flow quantitation weren't reported until the early 1980's (Moran et al., 1985; Nayler et al., 1986; Singer & Crooks, 1983; van Dijk, 1984). Current MR flow imaging methods are based on the fact that spins moving at a constant velocity accrue a phase proportional to the 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*<sup>0</sup> +*vt*. Rewriting equation 6, for *t* = *t*0:

$$
\phi = \gamma \int\_0^{t\_0} \vec{G}\_r(t) \cdot (\vec{r}\_0 + \vec{v}t) \, dt \tag{9}
$$

$$\vec{r} = \gamma \,\vec{r}\_0 \cdot \int\_0^{t\_0} \vec{G}\_r(t) \, dt + \gamma \,\vec{v} \cdot \int\_0^{t\_0} \vec{G}\_r(t) \, t \, dt \tag{10}$$

$$
\vec{r} = \gamma \vec{r}\_0 \cdot \vec{M}\_0 + \gamma \,\vec{v} \cdot \vec{M}\_1. \tag{11}
$$

where *M* <sup>0</sup> and *M* <sup>1</sup> are the zeroth and first moments of the*r*-gradient waveform at echo time, respectively. Thus, if a gradient with null zeroth moment is used (e.g., a bipolar gradient, aligned with *<sup>v</sup>*), the phase accrued for a constant velocity spin is *<sup>φ</sup>* <sup>=</sup> *<sup>γ</sup><sup>v</sup>* · *<sup>M</sup>* 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).

#### **2.3.1 Phase contrast**

6 Will-be-set-by-IN-TECH

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*(*kx*, *ky*) of the cross-sectional image

*m*(*x*, *y*)*e*

The Fourier coordinates *kx* and *ky* vary with time, according to the zeroth moment of the

 *t* 0

 *t* 0

These equations explain how the gradients can be used to "move" along k-space, as discussed on section 2.1. This formalism can be generalized for any combination of the gradients *Gx*, *Gy*

*m*(*r*) · *e*

 *t* 0

where *G <sup>r</sup>* is the oblique gradient resulting from the combination of the *Gx*, *Gy* and *Gz* gradients, and*r* is its corresponding axis, along which the linear variation in magnetic field

Given a spatial position function*r*(*t*) and a magnetic field gradient *G <sup>r</sup>*(*t*), the magnetization

 *t* 0

The basic principles of quantitative flow measurement using magnetic resonance were first proposed by Singer (1959) and Hahn (1960) in the late 1950's. However, clinical applications of MR flow quantitation weren't reported until the early 1980's (Moran et al., 1985; Nayler et al., 1986; Singer & Crooks, 1983; van Dijk, 1984). Current MR flow imaging methods are based on the fact that spins moving at a constant velocity accrue a phase proportional to the

 *t* 0

2*π*

2*π*

2*π*

<sup>−</sup>*j*2*π*(*kx <sup>x</sup>*+*ky <sup>y</sup>*) *dx dy*. (1)

*Gx*(*τ*) *dτ* (2)

*Gy*(*τ*) *dτ*. (3)

<sup>−</sup>*j*2*<sup>π</sup>kr*·*<sup>r</sup> <sup>d</sup><sup>r</sup>* (4)

*G <sup>r</sup>*(*τ*) *dτ*, (5)

*<sup>G</sup> <sup>r</sup>*(*τ*) ·*<sup>r</sup>*(*τ*) *<sup>d</sup>τ*, (6)

*G <sup>r</sup>*(*τ*) *dτ* (7)

<sup>=</sup> <sup>2</sup>*<sup>π</sup>kr* ·*<sup>r</sup>*, (8)

*M*(*kx*, *ky*) =

 *x y*

*kx*(*t*) = *<sup>γ</sup>*

*ky*(*t*) = *<sup>γ</sup>*

*kr*(*t*) = *<sup>γ</sup>*

*φ*(*r*, *t*) = *γ*

*φ* = *γr* ·

For static spins,*r*(*t*) is constant (*r*), and this becomes:

as in the exponential in equation 4.

**2.3 Principles of MR flow imaging**

*M*( *kr*) = *r*

**2.2 Mathematical formalism**

readout gradients *Gx* and *Gy*:

*m*(*x*, *y*):

and *Gz* as:

phase is:

intensity is perceived.

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:

$$w(\mathbf{x}, y) = \frac{\phi\_a(\mathbf{x}, y) - \phi\_b(\mathbf{x}, y)}{\gamma (M\_1^a - M\_1^b)}\,\tag{12}$$

where *φa*(*x*, *y*) and *φb*(*x*, *y*) are the phase images acquired in each acquisition, and *M<sup>a</sup>* <sup>1</sup> and *<sup>M</sup><sup>b</sup>* 1 are the first moment of the bipolar gradients used in each acquisition.

#### **2.3.2 Fourier velocity encoding**

While phase contrast provides a single velocity measurement associated with each voxel, Fourier velocity encoding (FVE) (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:

$$\phi(\vec{r}, \vec{v}, t) = 2\pi(\vec{k}\_r \cdot \vec{r} + \vec{k}\_\upsilon \cdot \vec{v}),\tag{13}$$

where*kv* is the velocity frequency variable associated with *v*, and is proportional to the first moment of *G <sup>r</sup>*(*t*):

$$
\vec{k}\_v = \frac{\gamma}{2\pi} \vec{M}\_1. \tag{14}
$$

**4.1 Pulse sequence**

*kv* encode level.

**4.2 Signal model**

approximately:

acquires one "disc" in *kx*-*ky*.

**RF**

**Gz**

**Gx**

**Gy**

timing corresponds to the studies shown in Figures 5 and 6.

kx

24:

1:

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 *kx*-*ky*-*kv* space (Figure 3). The bipolar gradient effectively phase-encodes in *kv*, while each spiral readout

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 11

**ab c d**

**1 ms 2.5 ms 8.1 ms 1.2 ms**

Fig. 2. Spiral FVE pulse sequence. It consists of (a) slice selective excitation, (b) velocity encoding bipolar gradient, (c) spiral readout, and (d) refocusing and spoiling gradients. This

kv

Fig. 3. Spiral FVE k-space sampling scheme. The dataset corresponding to each temporal frame is a stack-of-spirals in *kx*-*ky*-*kv* space. Each spiral acquisition corresponds to a different

2DFT phase contrast provides two two-dimensional functions, *m*(*x*, *y*) and *vo*(*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

ky

*s*(*x*, *y*, *v*) = *m*(*x*, *y*) × *δ*( *v* − *vo*(*x*, *y*) ), (15)

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, *M*(*kx*, *ky*, *kv*). Thus, an extra dimension is added to k-space, and multiple acquisitions are required to cover the entire *kx*-*ky*-*kv* space. In order to move along *kv*, a bipolar gradient with the appropriate amplitude (and first moment) is played before the *kx*-*ky* 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 *kv* 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 velocity distributions along the cross-sectional image *m*(*x*, *y*, *v*) is obtained by inverse Fourier transforming the acquired data *M*(*kx*, *ky*, *kv*). If cine imaging (Glover & Pelc, 1988) is used, measurements are also time resolved, resulting in a four-dimensional dataset: *m*(*x*, *y*, *v*, *t*).

The main drawback of FVE is scan time, as *kx*-*ky* should be fully sampled for each value of *kv*. As discussed in section 1.2.2, different approaches to accelerating FVE have been proposed. Those techniques are typically inefficient in spatially separating flowing blood from nearby static tissue. Furthermore, they are not capable of resolving multiple flows in a single acquisition. The spiral FVE method, proposed in section 4, addresses these limitations.

#### **3. Experimental setup**

Most experiments were performed on a Signa 3 T EXCITE HD system (GE Healthcare), with gradients capable of 40 mT/m amplitude and 150 T/m/s slew rate, and a receiver with sampling interval of 4 *μ*s. Sequence designs were optimized for this scanner configuration. The body coil was used for RF transmission in all studies. An 8-channel phase array cardiac coil was used in the healthy volunteer studies, but data from only 1 or 2 elements were used in reconstruction. In phantom studies, a single channel 5-inch surface coil was used. In the patient experiments presented in section 4, a Signa 1.5 T LX system (GE Healthcare) with the same gradient and receiver configuration was used, and acquisition was performed using a 5-inch surface coil.

The institutional review boards of the University of Southern California and Stanford University approved the imaging protocols. Subjects were screened for magnetic resonance imaging risk factors and provided informed consent in accordance with institutional policy.

#### **4. Slice-selective FVE with spiral readouts (spiral FVE)**

In order to address the limitations of existing flow imaging methods, 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 aortic valve, and comparisons with Doppler ultrasound and high-resolution 2DFT phase contrast MRI. The proposed method is demonstrated in healthy volunteers and in a patient.

#### **4.1 Pulse sequence**

8 Will-be-set-by-IN-TECH

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, *M*(*kx*, *ky*, *kv*). Thus, an extra dimension is added to k-space, and multiple acquisitions are required to cover the entire *kx*-*ky*-*kv* space. In order to move along *kv*, a bipolar gradient with the appropriate amplitude (and first moment) is played before the *kx*-*ky* 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 *kv* 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 velocity distributions along the cross-sectional image *m*(*x*, *y*, *v*) is obtained by inverse Fourier transforming the acquired data *M*(*kx*, *ky*, *kv*). If cine imaging (Glover & Pelc, 1988) is used, measurements are also time

The main drawback of FVE is scan time, as *kx*-*ky* should be fully sampled for each value of *kv*. As discussed in section 1.2.2, different approaches to accelerating FVE have been proposed. Those techniques are typically inefficient in spatially separating flowing blood from nearby static tissue. Furthermore, they are not capable of resolving multiple flows in a single acquisition. The spiral FVE method, proposed in section 4, addresses these limitations.

Most experiments were performed on a Signa 3 T EXCITE HD system (GE Healthcare), with gradients capable of 40 mT/m amplitude and 150 T/m/s slew rate, and a receiver with sampling interval of 4 *μ*s. Sequence designs were optimized for this scanner configuration. The body coil was used for RF transmission in all studies. An 8-channel phase array cardiac coil was used in the healthy volunteer studies, but data from only 1 or 2 elements were used in reconstruction. In phantom studies, a single channel 5-inch surface coil was used. In the patient experiments presented in section 4, a Signa 1.5 T LX system (GE Healthcare) with the same gradient and receiver configuration was used, and acquisition was performed using a

The institutional review boards of the University of Southern California and Stanford University approved the imaging protocols. Subjects were screened for magnetic resonance imaging risk factors and provided informed consent in accordance with institutional policy.

In order to address the limitations of existing flow imaging methods, 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.

We present practical implementations for measuring blood flow through the aortic valve, and comparisons with Doppler ultrasound and high-resolution 2DFT phase contrast MRI. The

resolved, resulting in a four-dimensional dataset: *m*(*x*, *y*, *v*, *t*).

**4. Slice-selective FVE with spiral readouts (spiral FVE)**

Scan-plane prescription is performed using classic protocols.

proposed method is demonstrated in healthy volunteers and in a patient.

**3. Experimental setup**

5-inch surface coil.

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 *kx*-*ky*-*kv* space (Figure 3). The bipolar gradient effectively phase-encodes in *kv*, while each spiral readout acquires one "disc" in *kx*-*ky*.

Fig. 2. Spiral FVE pulse sequence. It consists of (a) slice selective excitation, (b) velocity encoding bipolar gradient, (c) spiral readout, and (d) refocusing and spoiling gradients. This timing corresponds to the studies shown in Figures 5 and 6.

Fig. 3. Spiral FVE k-space sampling scheme. The dataset corresponding to each temporal frame is a stack-of-spirals in *kx*-*ky*-*kv* space. Each spiral acquisition corresponds to a different *kv* encode level.

#### **4.2 Signal model**

2DFT phase contrast provides two two-dimensional functions, *m*(*x*, *y*) and *vo*(*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:

$$s(\mathbf{x}, y, \upsilon) = m(\mathbf{x}, y) \times \delta(\upsilon - v\_o(\mathbf{x}, y) \text{ )},\tag{15}$$

**healthy volunteer patient**

**field-of-view** 25 cm 20 cm **spatial resolution** 7 mm 6.5 mm

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 13

*kv* **encodes** 24 64 **velocity field-of-view** 600 or 800 cm/s 1200 cm/s **velocity resolution** 25 or 33 cm/s 19 cms/s

*kv* **encodes/heartbeat** 2 4 **TR** 12.8 or 12.5 ms 12.8 ms **temporal resolution** 25.6 or 25 ms 51.2 ms

Table 1. Scan parameters used in the different spiral FVE studies.

using the same data (see Figure 8, discussed later).

40 minutes for phase contrast, and 12 seconds for FVE.

**4.5 Accuracy of spiral FVE measurements**

high-resolution 2DFT phase contrast data.

**scan time** 12 heartbeats 16 heartbeats

then applied to *S*ROI(*kv*, *t*) to increase the number of temporal frames (Riederer et al., 1988). Saturation effects (Gao et al., 1988) are compensated by normalizing the 2-norm of *S*ROI(*kv*) independently for each temporal frame, which effectively normalizes each cardiac phase. *S*ROI(*kv*, *t*) is then zero-padded along the *kv* axis, and an inverse Fourier transform produces *s*ROI(*v*, *t*). The time-velocity histogram for the ROI is |*s*ROI(*v*, *t*)|, and for display purposes, smoother histograms are obtained by cubic spline interpolation along *t* (Bartels et al., 1987). The reconstruction process can be repeated for each voxel, or for multiple regions of interest,

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 on section 4.2 was used to generate simulated FVE data based on

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* = 3 mm and Δ*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

A simulated spiral FVE dataset was computed from the PC magnitude and velocity maps, using the convolution model described in Eq. 16. 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*(*kx*, *ky*) 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.

where *δ*(*v*) is the Dirac delta function.

As spiral FVE acquisitions follow a stack-of-spirals pattern in *kx*-*ky*-*kv* space (Figure 3), k-space data is truncated to a cylinder, i.e., a circle along *kx*-*ky* (with diameter 1/Δ*r*), and a rectangle along *kv* (with width 1/Δ*v*), where Δ*r* and Δ*v* are the prescribed spatial and velocity resolutions, respectively. The associated object domain spatial-velocity blurring can be modeled as a convolution of the true object distribution, *s*(*x*, *y*, *v*), with jinc( *x*<sup>2</sup> + *y*2/Δ*r*) and sinc(*v*/Δ*v*), resulting in:

$$\hat{s}(\mathbf{x}, \boldsymbol{y}, \upsilon) = \left[ m(\mathbf{x}, \boldsymbol{y}) \times \delta(\upsilon - \upsilon\_o(\mathbf{x}, \boldsymbol{y}) \,) \right] \ast \text{sinc} \left( \frac{\upsilon}{\Delta \upsilon} \right) \ast \text{jinc} \left( \frac{\sqrt{\mathbf{x}^2 + \boldsymbol{y}^2}}{\Delta r} \right)$$

$$= \left[ m(\mathbf{x}, \boldsymbol{y}) \times \text{sinc} \left( \frac{\upsilon - \upsilon\_o(\mathbf{x}, \boldsymbol{y})}{\Delta \upsilon} \right) \right] \ast \text{jinc} \left( \frac{\sqrt{\mathbf{x}^2 + \boldsymbol{y}^2}}{\Delta r} \right) \,. \tag{16}$$

where *s*ˆ(*x*, *y*, *v*) is the measured object distribution, and ∗ denotes convolution.

#### **4.3 Data acquisition**

The excitation achieved a 5 mm slice thickness and 30◦ flip angle, with a 0.5 ms RF pulse and a 1 ms gradient. Through-plane velocity encoding was implemented using a large bipolar pulse along the *z* direction that was scaled to achieve different *kv* encodes. The velocity resolution is determined by the first moment of the largest bipolar gradient, and the velocity field-of-view is determined by the increment in gradient first moment for different velocity encodes. A bipolar duration of 2.5 ms achieves a velocity resolution of 25 cm/s. Gradient duration increases if velocity resolution is improved.

An 8.1 ms optimized uniform-density single-shot spiral acquisition (Hargreaves, 2001) acquires *kx*-*ky* at each velocity encode, and zeroth and first moments are refocused in 0.5 ms. The readout and refocusing gradients were designed using public domain software1. A 0.65 ms spoiling gradient (Zur et al., 1991) achieves a 6*π* phase-wrap over the slice thickness. The spoiling gradient was not overlapped with the refocusing gradients, but this could be done to further shorten the TR. The minimum TR (approximately 13 ms) was used in all studies. Other scan dependent pulse sequence parameters are listed in Table 1.

Prospective ECG gating was used to synchronize acquisitions with the cardiac cycle. Two *kv* levels were repeatedly acquired during each heartbeat in order to resolve 25 to 35 cardiac phases and produce a cine dataset (Figure 5, discussed later). The true temporal resolution was 26 ms (2 TRs). Sliding window reconstruction (Riederer et al., 1988) was used to produce a new image every 13 ms.

#### **4.4 Data reconstruction**

Reconstruction was performed in Matlab (The MathWorks, Inc., Natick, MA, USA). Each spiral interleaf is first gridded (Jackson et al., 1991) and inverse Fourier transformed to form an image, *m*(*x*, *y*), for each temporal frame. This step converts the acquired data *Skx*,*ky* (*kv*, *t*) to *Sx*,*y*(*kv*, *t*). The operator manually defines a region of interest (ROI) in the *x*-*y* plane using the image corresponding to *kv* = 0 and *t* = 0. Pixel intensities within the ROI are averaged at each temporal frame, resulting in a 2D dataset: *S*ROI(*kv*, *t*) = ∑ROI *<sup>x</sup>*,*<sup>y</sup> Sx*,*y*(*kv*, *t*). View sharing is

<sup>1</sup> http://www-mrsrl.stanford.edu/∼brian/vdspiral/



Table 1. Scan parameters used in the different spiral FVE studies.

then applied to *S*ROI(*kv*, *t*) to increase the number of temporal frames (Riederer et al., 1988). Saturation effects (Gao et al., 1988) are compensated by normalizing the 2-norm of *S*ROI(*kv*) independently for each temporal frame, which effectively normalizes each cardiac phase. *S*ROI(*kv*, *t*) is then zero-padded along the *kv* axis, and an inverse Fourier transform produces *s*ROI(*v*, *t*). The time-velocity histogram for the ROI is |*s*ROI(*v*, *t*)|, and for display purposes, smoother histograms are obtained by cubic spline interpolation along *t* (Bartels et al., 1987). The reconstruction process can be repeated for each voxel, or for multiple regions of interest, using the same data (see Figure 8, discussed later).

#### **4.5 Accuracy of spiral FVE measurements**

As spiral FVE acquisitions follow a stack-of-spirals pattern in *kx*-*ky*-*kv* space (Figure 3), k-space data is truncated to a cylinder, i.e., a circle along *kx*-*ky* (with diameter 1/Δ*r*), and a rectangle along *kv* (with width 1/Δ*v*), where Δ*r* and Δ*v* are the prescribed spatial and velocity resolutions, respectively. The associated object domain spatial-velocity blurring can

> *<sup>v</sup>* <sup>−</sup> *vo*(*x*, *<sup>y</sup>*) Δ*v*

The excitation achieved a 5 mm slice thickness and 30◦ flip angle, with a 0.5 ms RF pulse and a 1 ms gradient. Through-plane velocity encoding was implemented using a large bipolar pulse along the *z* direction that was scaled to achieve different *kv* encodes. The velocity resolution is determined by the first moment of the largest bipolar gradient, and the velocity field-of-view is determined by the increment in gradient first moment for different velocity encodes. A bipolar duration of 2.5 ms achieves a velocity resolution of 25 cm/s. Gradient

An 8.1 ms optimized uniform-density single-shot spiral acquisition (Hargreaves, 2001) acquires *kx*-*ky* at each velocity encode, and zeroth and first moments are refocused in 0.5 ms. The readout and refocusing gradients were designed using public domain software1. A 0.65 ms spoiling gradient (Zur et al., 1991) achieves a 6*π* phase-wrap over the slice thickness. The spoiling gradient was not overlapped with the refocusing gradients, but this could be done to further shorten the TR. The minimum TR (approximately 13 ms) was used in all

Prospective ECG gating was used to synchronize acquisitions with the cardiac cycle. Two *kv* levels were repeatedly acquired during each heartbeat in order to resolve 25 to 35 cardiac phases and produce a cine dataset (Figure 5, discussed later). The true temporal resolution was 26 ms (2 TRs). Sliding window reconstruction (Riederer et al., 1988) was used to produce

Reconstruction was performed in Matlab (The MathWorks, Inc., Natick, MA, USA). Each spiral interleaf is first gridded (Jackson et al., 1991) and inverse Fourier transformed to form an image, *m*(*x*, *y*), for each temporal frame. This step converts the acquired data *Skx*,*ky* (*kv*, *t*) to *Sx*,*y*(*kv*, *t*). The operator manually defines a region of interest (ROI) in the *x*-*y* plane using the image corresponding to *kv* = 0 and *t* = 0. Pixel intensities within the ROI are averaged at

 *v* Δ*v* ∗ jinc

∗ jinc

*x*<sup>2</sup> + *y*2/Δ*r*)

*<sup>x</sup>*,*<sup>y</sup> Sx*,*y*(*kv*, *t*). View sharing is

, (16)

*<sup>x</sup>*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*<sup>2</sup> Δ*r*

*<sup>x</sup>*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*<sup>2</sup> Δ*r*

be modeled as a convolution of the true object distribution, *s*(*x*, *y*, *v*), with jinc(

where *s*ˆ(*x*, *y*, *v*) is the measured object distribution, and ∗ denotes convolution.

studies. Other scan dependent pulse sequence parameters are listed in Table 1.

each temporal frame, resulting in a 2D dataset: *S*ROI(*kv*, *t*) = ∑ROI

<sup>1</sup> http://www-mrsrl.stanford.edu/∼brian/vdspiral/

*s*ˆ(*x*, *y*, *v*) = [*m*(*x*, *y*) × *δ*( *v* − *vo*(*x*, *y*) )] ∗ sinc

*m*(*x*, *y*) × sinc

duration increases if velocity resolution is improved.

where *δ*(*v*) is the Dirac delta function.

= 

and sinc(*v*/Δ*v*), resulting in:

**4.3 Data acquisition**

a new image every 13 ms.

**4.4 Data reconstruction**

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 on section 4.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* = 3 mm and Δ*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. 16. 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*(*kx*, *ky*) 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.

• • • 26 ms

trigger

HB #1 1212 1212 • • •

HB #2 3434 3434 • • • HB #3 5656 5656 • • •

HB #12 23 24 23 24 23 24 23 24 • • •

levels are acquired in this sequential fashion (see white arrow).

0 200 400 600 800

**FVE** time (ms)

Fig. 6. Comparison of the spiral FVE method with Doppler ultrasound, in a healthy

time (ms)

**Ultrasound**

**Ultrasound**

0 200 400 600 800 −300

−150

0

150

velocity (cm/s)

Fig. 5. Artifacts with sequential view-ordering scheme. Each box represents the acquisition of one *kv* level, during one imaging TR. A sliding window is used to produce a new image every TR. Ghosting artifacts appear shifted by 1/2 of the velocity field-of-view when the *kv*

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 15

Artifacts and loss of temporal resolution due to view sharing can be avoided or corrected using different approaches. Acquiring multiple *kv* levels per heartbeat reduces scan time, but causes blurring along the time axis and ghosting along the velocity axis. Blurring is caused by the reduction in temporal resolution, and ghosting artifacts arise when the velocity distribution changes between the acquisition of consecutive velocity encodes. Both ghosting and blurring can be overcome by acquiring only one view per heartbeat, but this would require increase in scan time or reduction in velocity resolution. As an alternative, these artifacts may be corrected using techniques that exploit efficient use of *k*-*t* space, such as the

A representative *in vivo* result is compared with Doppler ultrasound in Figure 6. The MRI measured time-velocity histogram show good agreement with the ultrasound measurement, as the peak velocity and the shape of the flow waveform were comparable to those observed

**FVE**

c

300

Frame #2 Frame #4

Frame #1 Frame #3

approach proposed in section 5.

in the ultrasound study.

velocity (cm/s)

−50

volunteer aortic valve study.

50

150

**4.6.2 Spiral FVE vs. doppler ultrasound**

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 measured with spiral FVE agree well with those obtained with 2DFT phase contrast, the current MRI gold standard. This approach for deriving FVE data from high-resolution velocity maps (Eq. 16) can be used for many simulation purposes (Carvalho et al., 2010).

10.7 dB 11.3 dB 9.3 dB 10.8 dB 10.3 dB 10.5 dB 10.8 dB 10.9 dB 11.7 dB **e**

Fig. 4. *In vitro* evaluation of the accuracy of spiral FVE velocity histograms. Results are shown for nine representative voxels, selected around the circumference of the vessel wall of the pulsatile carotid flow phantom's bifurcation (a). For each voxel, it is shown: (b) time-velocity distribution derived from high-resolution 2DFT phase contrast; (c) time-velocity distribution measured with spiral FVE; (d) absolute difference between spiral FVE and 2DFT PC-derived histograms; (e) signal-to-error ratio.

#### **4.6 Aortic valve flow assessment using spiral FVE**

The proposed method was evaluated *in vivo*, aiming at quantifying flow through the aortic valve. Scan-plane prescription was performed using a real-time imaging sequence.

For a severely stenosed heart valve, peak velocities can reach up to 600 cm/s (Galea et al., 2002). As regurgitant jets don't overlap in time with forward flow, we used a 600 cm/s velocity field-of-view (±300 cm/s). This value could be increased by extending the scan time, or by sacrificing temporal, velocity and/or spatial resolutions (see Figure 9, discussed later). Scan parameters were summarized in Table 1.

#### **4.6.1 Order of velocity encode acquisitions**

When the *kv* levels are acquired in a sequential fashion, ghosting artifacts due to data inconsistency appear shifted by one half of the velocity field-of-view (Figure 5). Using this sampling scheme and an appropriate velocity field-of-view, the artifacts will not overlap with the flow profile and may be easily identified and masked out.

Fig. 5. Artifacts with sequential view-ordering scheme. Each box represents the acquisition of one *kv* level, during one imaging TR. A sliding window is used to produce a new image every TR. Ghosting artifacts appear shifted by 1/2 of the velocity field-of-view when the *kv* levels are acquired in this sequential fashion (see white arrow).

Artifacts and loss of temporal resolution due to view sharing can be avoided or corrected using different approaches. Acquiring multiple *kv* levels per heartbeat reduces scan time, but causes blurring along the time axis and ghosting along the velocity axis. Blurring is caused by the reduction in temporal resolution, and ghosting artifacts arise when the velocity distribution changes between the acquisition of consecutive velocity encodes. Both ghosting and blurring can be overcome by acquiring only one view per heartbeat, but this would require increase in scan time or reduction in velocity resolution. As an alternative, these artifacts may be corrected using techniques that exploit efficient use of *k*-*t* space, such as the approach proposed in section 5.

#### **4.6.2 Spiral FVE vs. doppler ultrasound**

12 Will-be-set-by-IN-TECH

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 measured with spiral FVE agree well with those obtained with 2DFT phase contrast, the current MRI gold standard. This approach for deriving FVE data from high-resolution velocity maps (Eq. 16) can be used for many simulation purposes (Carvalho et al., 2010).

> **b** velocity

1

**c**

**d**

**e**

**4.6 Aortic valve flow assessment using spiral FVE**

parameters were summarized in Table 1.

**4.6.1 Order of velocity encode acquisitions**

10.7 dB

FVE and 2DFT PC-derived histograms; (e) signal-to-error ratio.

the flow profile and may be easily identified and masked out.

11.3 dB

9.3 dB

Fig. 4. *In vitro* evaluation of the accuracy of spiral FVE velocity histograms. Results are shown for nine representative voxels, selected around the circumference of the vessel wall of

time-velocity distribution measured with spiral FVE; (d) absolute difference between spiral

The proposed method was evaluated *in vivo*, aiming at quantifying flow through the aortic

For a severely stenosed heart valve, peak velocities can reach up to 600 cm/s (Galea et al., 2002). As regurgitant jets don't overlap in time with forward flow, we used a 600 cm/s velocity field-of-view (±300 cm/s). This value could be increased by extending the scan time, or by sacrificing temporal, velocity and/or spatial resolutions (see Figure 9, discussed later). Scan

When the *kv* levels are acquired in a sequential fashion, ghosting artifacts due to data inconsistency appear shifted by one half of the velocity field-of-view (Figure 5). Using this sampling scheme and an appropriate velocity field-of-view, the artifacts will not overlap with

the pulsatile carotid flow phantom's bifurcation (a). For each voxel, it is shown: (b) time-velocity distribution derived from high-resolution 2DFT phase contrast; (c)

valve. Scan-plane prescription was performed using a real-time imaging sequence.

10.8 dB

10.3 dB

10.5 dB

10.8 dB

10.9 dB

11.7 dB

1

**a**

2 3 4

9 8 7

6 5

time

2

3

4

5

6

7

8

9

A representative *in vivo* result is compared with Doppler ultrasound in Figure 6. The MRI measured time-velocity histogram show 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 study.

**Ultrasound**

Fig. 6. Comparison of the spiral FVE method with Doppler ultrasound, in a healthy volunteer aortic valve study.

resolution.

the breath-hold duration to acquire more *kv* levels, or by reducing the velocity field-of-view. Temporal resolution can be made as high as one TR duration (13 ms) by segmenting the *kv* encodes across additional heartbeats (longer breath-holds), or by compromising velocity resolution or field-of-view. Spatial resolution can be improved by reducing the spatial field-of-view, or by increasing the number of spiral interleaves, which would require compromising other scan parameters such as scan time, temporal resolution and/or velocity

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 17

10 30 50 70

temporal resolution (ms)

breath-hold duration. Velocity resolution corresponds to a 600 cm/s field-of-view, temporal resolution corresponds to a 8.1 ms spiral readout, and scan time corresponds to a single-shot spiral acquisition. The arrow indicates the configuration used in the study in Figure 6 (2·TR

As the spiral readouts are considerably long, a potential issue in spiral FVE imaging is blurring in image domain, due to off-resonance. Because SNR was not a limiting issue for the applications we have presented, spiral FVE may perform better at lower field strengths where there is reduced off-resonance. At 3 T, localized shimming and off-resonance correction techniques can be used to reduce blurring. Furthermore, readout duration can be reduced by decreasing the spatial resolution or field-of-view, or by using variable-density spirals (Tsai & Nishimura, 2000). Another alternative is to use multiple short spiral interleaves, which would require longer scan times, but parallel imaging techniques (Pruessmann et al., 2001; Samsonov et al., 2006) can potentially accelerate acquisition if multi-channel receiver coils are used. This approach also has the benefit of allowing increase in the frame rate, as the number of imaged cardiac phases is limited by the minimum TR. Another possible solution to the off-resonance problem is the use of echo-planar imaging trajectories, which produce different off-resonance effects (geometric warping) (Feinberg & Oshio, 1992), but are also more sensitive to artifacts

A noticeable artifact in spiral FVE is Gibbs ringing along the velocity dimension. These artifacts can be less noticeable if velocity resolution is increased, which would also improve the ability to visualize features in the flow waveform and the precision to resolve the peak velocity, but would require longer breath-holds. Alternatively, the velocity resolution can be improved by using variable-density sampling along *kv* (Carvalho et al., 2007; DiCarlo et al., 2005), or partial Fourier techniques (Noll, Nishimura & Macovski, 1991). Another approach to reducing ringing artifacts is to window the *kv* samples before applying the inverse

Fig. 9. Spiral FVE trade-offs between temporal resolution, velocity resolution, and

8 heartbeats 12 heartbeats 16 heartbeats 24 heartbeats 32 heartbeats

0

20

40

velocity resolution (cm/s)

temporal resolution, 24 velocity encodes).

**4.8 Issues with spiral FVE**

from in-plane flow or motion.

60

80

#### **4.6.3 Patient evaluation**

Figure 7 shows the time-velocity distribution measured through the aortic valve of a patient with aortic stenosis. This result demonstrate that spiral FVE can accurately detect complex flow, as a high-speed jet with a wide distribution of velocities is clearly visible.

Fig. 7. Evaluation of spiral FVE in a patient with aortic stenosis. Note the high-speed jet with a wide distribution of velocities.

#### **4.6.4 Measurement of multiple flows**

Figure 8 illustrates spiral FVE's ability of resolving different flows from a single dataset. A different flow distribution was calculated for each voxel, and the distributions from single voxels from different ROIs are displayed. Red and blue dots indicate voxels where ascending and descending blood flows were detected, respectively, and the color intensity of each dot indicates the highest velocity detected in that voxel in a particular temporal frame.

Fig. 8. Multiple flow distributions obtained from a single spiral FVE dataset. For each voxel in the image, a flow distribution was calculated, and the red and blue dots indicate voxels where ascending and descending blood flows were detected, respectively. The color intensity of each dot indicates the highest velocity detected in that voxel in a particular temporal frame (indicated by the white dashed lines).

#### **4.7 Resolution trade-offs**

In spiral FVE, there is an important trade-off between velocity resolution, temporal resolution, and scan time (Figure 9). This trade-off also involves other scan parameters, such as velocity field-of-view, number of spiral interleaves, spiral readout duration, spatial resolution, and spatial field-of-view. Velocity resolution can be improved in many ways, such as increasing the breath-hold duration to acquire more *kv* levels, or by reducing the velocity field-of-view. Temporal resolution can be made as high as one TR duration (13 ms) by segmenting the *kv* encodes across additional heartbeats (longer breath-holds), or by compromising velocity resolution or field-of-view. Spatial resolution can be improved by reducing the spatial field-of-view, or by increasing the number of spiral interleaves, which would require compromising other scan parameters such as scan time, temporal resolution and/or velocity resolution.

Fig. 9. Spiral FVE trade-offs between temporal resolution, velocity resolution, and breath-hold duration. Velocity resolution corresponds to a 600 cm/s field-of-view, temporal resolution corresponds to a 8.1 ms spiral readout, and scan time corresponds to a single-shot spiral acquisition. The arrow indicates the configuration used in the study in Figure 6 (2·TR temporal resolution, 24 velocity encodes).

#### **4.8 Issues with spiral FVE**

14 Will-be-set-by-IN-TECH

Figure 7 shows the time-velocity distribution measured through the aortic valve of a patient with aortic stenosis. This result demonstrate that spiral FVE can accurately detect complex

time (ms)

0 100 200 300 400 500 −200

Fig. 7. Evaluation of spiral FVE in a patient with aortic stenosis. Note the high-speed jet with

Figure 8 illustrates spiral FVE's ability of resolving different flows from a single dataset. A different flow distribution was calculated for each voxel, and the distributions from single voxels from different ROIs are displayed. Red and blue dots indicate voxels where ascending and descending blood flows were detected, respectively, and the color intensity of each dot

Fig. 8. Multiple flow distributions obtained from a single spiral FVE dataset. For each voxel in the image, a flow distribution was calculated, and the red and blue dots indicate voxels where ascending and descending blood flows were detected, respectively. The color intensity of each dot indicates the highest velocity detected in that voxel in a particular temporal

In spiral FVE, there is an important trade-off between velocity resolution, temporal resolution, and scan time (Figure 9). This trade-off also involves other scan parameters, such as velocity field-of-view, number of spiral interleaves, spiral readout duration, spatial resolution, and spatial field-of-view. Velocity resolution can be improved in many ways, such as increasing

<sup>0</sup> <sup>400</sup> <sup>800</sup> −200

0

−150

0

150

200

indicates the highest velocity detected in that voxel in a particular temporal frame.

flow, as a high-speed jet with a wide distribution of velocities is clearly visible.

velocity (cm/s)

0

200

400

**4.6.3 Patient evaluation**

a wide distribution of velocities.

**4.6.4 Measurement of multiple flows**

<sup>0</sup> <sup>400</sup> <sup>800</sup> −200

<sup>0</sup> <sup>400</sup> <sup>800</sup> −200

frame (indicated by the white dashed lines).

0

0

**4.7 Resolution trade-offs**

200

200

As the spiral readouts are considerably long, a potential issue in spiral FVE imaging is blurring in image domain, due to off-resonance. Because SNR was not a limiting issue for the applications we have presented, spiral FVE may perform better at lower field strengths where there is reduced off-resonance. At 3 T, localized shimming and off-resonance correction techniques can be used to reduce blurring. Furthermore, readout duration can be reduced by decreasing the spatial resolution or field-of-view, or by using variable-density spirals (Tsai & Nishimura, 2000). Another alternative is to use multiple short spiral interleaves, which would require longer scan times, but parallel imaging techniques (Pruessmann et al., 2001; Samsonov et al., 2006) can potentially accelerate acquisition if multi-channel receiver coils are used. This approach also has the benefit of allowing increase in the frame rate, as the number of imaged cardiac phases is limited by the minimum TR. Another possible solution to the off-resonance problem is the use of echo-planar imaging trajectories, which produce different off-resonance effects (geometric warping) (Feinberg & Oshio, 1992), but are also more sensitive to artifacts from in-plane flow or motion.

A noticeable artifact in spiral FVE is Gibbs ringing along the velocity dimension. These artifacts can be less noticeable if velocity resolution is increased, which would also improve the ability to visualize features in the flow waveform and the precision to resolve the peak velocity, but would require longer breath-holds. Alternatively, the velocity resolution can be improved by using variable-density sampling along *kv* (Carvalho et al., 2007; DiCarlo et al., 2005), or partial Fourier techniques (Noll, Nishimura & Macovski, 1991). Another approach to reducing ringing artifacts is to window the *kv* samples before applying the inverse

**5.1.2 Acceleration via UNFOLD**

notch filter along *t*.

**a b**

**0**

**kv (ms/cm)**

**14.2**

**−15**

**interleaf 1 interleaf 2 interleaf 3**

reduces overlaps and facilitates filtering.

Scan time in FVE imaging can be significantly reduced using multi-dimensional temporal acceleration (Gamper et al., 2008; Hansen et al., 2004). An implementation of the UNFOLD method (Madore et al., 1999; Tsao, 2002) was specially designed for spiral FVE. A view-ordering scheme that reduces overlap in *v*-*f* space (*f* denotes temporal frequency) was designed. It consists in alternating spiral interleaves and *kv* encodings for each cardiac phase, according to Figure 10a. The associated point spread function is such that aliasing replicas, caused by temporal undersampling, are separated from the main lobe both in velocity (by half of the velocity FOV) and in temporal frequency (by 1/2TR) (Figure 10b) (Hansen et al., 2004; Tsao, 2002). Aliasing components caused by the sidelobes at ±20 and ±40 Hz are expected to correspond to static or slow moving spins, and hence will have a small footprint in *v*-*f* space. This is because these sidelobes spread around the periphery of the spatial FOV, but are null at the center, where high-velocity pulsatile flow is located. The aliasing signal is filtered using the two-dimensional filter shown in Figure 11. This filter has a bandwidth of 107 Hz for velocities below ±150 cm/s. For higher velocities, the bandwidth varies from 69 to 30 Hz. This results in effective temporal resolutions of 9.3 ms and 14.5–33.3 ms, respectively. The temporal resolution is lower for higher velocities, but this may prove unnoticeable, as the velocity distribution of high-velocity flow jets within large voxels is typically temporally smooth. For comparison, the temporal resolution with the conventional approach — view sharing (Riederer et al., 1988) — would be 50 ms for all *v*. The remaining narrow-bandwidth aliasing components at ±20 and ±40 Hz are filtered using a tight zero-phase one-dimensional

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 19

**time**

400

600

200

0

velocity (cm/s)

−200

−400

−600

Fig. 10. Proposed view-ordering scheme for accelerated spiral FVE (a) and its corresponding point spread function (b). In (a), each color represents a different spiral interleaf, and darker tones indicate "views" that are discarded in the partial Fourier experiments. Views aligned in *kv* are acquired sequentially throughout the cardiac cycle. Views aligned in time (same cardiac phase) are acquired in different heartbeats. In (b), each square shows the point spread function in *x*-*y* for a particular *v*-*f* coordinate. Aliasing replicas are separated from the main lobe by half of the velocity FOV, and half of the temporal frequency bandwidth, which

−60

−40 −20 <sup>0</sup> <sup>20</sup> <sup>40</sup> −40 dB

25 cm

25 cm

temporal frequency (Hz)

−35 dB −30 dB −25 dB −20 dB −15 dB −10 dB −5 dB 0 dB

60

y x

Fourier transform (Bernstein et al., 2001). However, windows with lower sidelobes generally have wider mainlobes, which would cause blurring along the velocity axis and consequent reduction in velocity resolution.

One drawback of the proposed method is the requirement of cardiac gating and breath-holding. Cardiac gating does not work well in patients with arrhythmias, and breath-holding may cause hemodynamic changes and is not possible for some patients (Macgowan et al., 2005). However, arrhythmia rejection (Chia et al., 2000) and respiratory gating schemes may overcome these problems, at the cost of increased scan time.

## **5. Accelerated spiral FVE**

As introduced in section 4, spiral FVE presents limitations such as insufficient velocity field-of-view (FOV), low spatial resolution, and moderate temporal resolution. In particular, the use of view sharing (Riederer et al., 1988) causes blurring along the temporal dimension (*t*) and ghosting along the velocity dimension (*v*) (Figure 5), and the use of long spiral readouts makes the technique sensitive to off-resonance, resulting in blurring along the in-plane spatial axes (*x* and *y*). The approach proposed in this section aims to address these limitations.

Spiral FVE datasets are four-dimensional, which makes this method particularly suitable for accelerated acquisition (Hansen et al., 2004). In this section, we achieve 18-fold acceleration using a combination of three techniques: variable-density spiral sampling along *kx*-*ky* (Tsai & Nishimura, 2000); partial Fourier (Noll, Nishimura & Macovski, 1991) along *kv*; and temporal acceleration through a novel implementation of the UNFOLD method (Madore et al., 1999; Tsao, 2002). The improved acquisition is performed without increase in scan time compared with the original implementation, and is demonstrated *in vivo* in a healthy volunteer.

#### **5.1 Accelerated data acquisition**

#### **5.1.1 Acceleration via variable-density sampling**

Variable-density spirals have been shown to increase spatiotemporal resolution and improve accuracy in flow quantitation (Liu et al., 2008). The spatial aliasing resulting from variable-density spiral sampling is incoherent, and, in the regions-of-interest (e.g., cardiac chambers, valves, great vessels), it typically originates from static or slow moving material located at the periphery of the spatial FOV (e.g., chest wall). FVE resolves the distribution of velocities within the voxel, thus moderate low-velocity aliasing artifacts generally do not affect one's ability to calculate diagnostically important parameters — such as peak velocity and acceleration — from the time-velocity distribution. Spiral FVE's single-shot uniform-density spiral readout was replaced with a multi-shot variable-density spiral acquisition (Tsai & Nishimura, 2000). The use of multi-shot acquisitions provides the possibility of multi-dimensional temporal acceleration, and allows reduction of readout duration and TR, which reduce off-resonance artifacts and temporal aliasing, respectively. The use of a shorter TR also allows improving the temporal resolution. Gradient waveforms were designed using public domain software2, based on the hardware limits of our scanner. The spatial FOV was varied linearly from 25 cm at the center of k-space to 6.25 cm at the periphery.

<sup>2</sup> http://www-mrsrl.stanford.edu/∼brian/vdspiral/

#### **5.1.2 Acceleration via UNFOLD**

16 Will-be-set-by-IN-TECH

Fourier transform (Bernstein et al., 2001). However, windows with lower sidelobes generally have wider mainlobes, which would cause blurring along the velocity axis and consequent

One drawback of the proposed method is the requirement of cardiac gating and breath-holding. Cardiac gating does not work well in patients with arrhythmias, and breath-holding may cause hemodynamic changes and is not possible for some patients (Macgowan et al., 2005). However, arrhythmia rejection (Chia et al., 2000) and respiratory gating schemes may overcome these problems, at the cost of increased scan time.

As introduced in section 4, spiral FVE presents limitations such as insufficient velocity field-of-view (FOV), low spatial resolution, and moderate temporal resolution. In particular, the use of view sharing (Riederer et al., 1988) causes blurring along the temporal dimension (*t*) and ghosting along the velocity dimension (*v*) (Figure 5), and the use of long spiral readouts makes the technique sensitive to off-resonance, resulting in blurring along the in-plane spatial axes (*x* and *y*). The approach proposed in this section aims to address these limitations. Spiral FVE datasets are four-dimensional, which makes this method particularly suitable for accelerated acquisition (Hansen et al., 2004). In this section, we achieve 18-fold acceleration using a combination of three techniques: variable-density spiral sampling along *kx*-*ky* (Tsai & Nishimura, 2000); partial Fourier (Noll, Nishimura & Macovski, 1991) along *kv*; and temporal acceleration through a novel implementation of the UNFOLD method (Madore et al., 1999; Tsao, 2002). The improved acquisition is performed without increase in scan time compared

with the original implementation, and is demonstrated *in vivo* in a healthy volunteer.

Variable-density spirals have been shown to increase spatiotemporal resolution and improve accuracy in flow quantitation (Liu et al., 2008). The spatial aliasing resulting from variable-density spiral sampling is incoherent, and, in the regions-of-interest (e.g., cardiac chambers, valves, great vessels), it typically originates from static or slow moving material located at the periphery of the spatial FOV (e.g., chest wall). FVE resolves the distribution of velocities within the voxel, thus moderate low-velocity aliasing artifacts generally do not affect one's ability to calculate diagnostically important parameters — such as peak velocity and acceleration — from the time-velocity distribution. Spiral FVE's single-shot uniform-density spiral readout was replaced with a multi-shot variable-density spiral acquisition (Tsai & Nishimura, 2000). The use of multi-shot acquisitions provides the possibility of multi-dimensional temporal acceleration, and allows reduction of readout duration and TR, which reduce off-resonance artifacts and temporal aliasing, respectively. The use of a shorter TR also allows improving the temporal resolution. Gradient waveforms were designed using public domain software2, based on the hardware limits of our scanner. The spatial FOV was varied linearly from 25 cm at the center of k-space to 6.25 cm at the

reduction in velocity resolution.

**5. Accelerated spiral FVE**

**5.1 Accelerated data acquisition**

periphery.

**5.1.1 Acceleration via variable-density sampling**

<sup>2</sup> http://www-mrsrl.stanford.edu/∼brian/vdspiral/

Scan time in FVE imaging can be significantly reduced using multi-dimensional temporal acceleration (Gamper et al., 2008; Hansen et al., 2004). An implementation of the UNFOLD method (Madore et al., 1999; Tsao, 2002) was specially designed for spiral FVE. A view-ordering scheme that reduces overlap in *v*-*f* space (*f* denotes temporal frequency) was designed. It consists in alternating spiral interleaves and *kv* encodings for each cardiac phase, according to Figure 10a. The associated point spread function is such that aliasing replicas, caused by temporal undersampling, are separated from the main lobe both in velocity (by half of the velocity FOV) and in temporal frequency (by 1/2TR) (Figure 10b) (Hansen et al., 2004; Tsao, 2002). Aliasing components caused by the sidelobes at ±20 and ±40 Hz are expected to correspond to static or slow moving spins, and hence will have a small footprint in *v*-*f* space. This is because these sidelobes spread around the periphery of the spatial FOV, but are null at the center, where high-velocity pulsatile flow is located. The aliasing signal is filtered using the two-dimensional filter shown in Figure 11. This filter has a bandwidth of 107 Hz for velocities below ±150 cm/s. For higher velocities, the bandwidth varies from 69 to 30 Hz. This results in effective temporal resolutions of 9.3 ms and 14.5–33.3 ms, respectively. The temporal resolution is lower for higher velocities, but this may prove unnoticeable, as the velocity distribution of high-velocity flow jets within large voxels is typically temporally smooth. For comparison, the temporal resolution with the conventional approach — view sharing (Riederer et al., 1988) — would be 50 ms for all *v*. The remaining narrow-bandwidth aliasing components at ±20 and ±40 Hz are filtered using a tight zero-phase one-dimensional notch filter along *t*.

Fig. 10. Proposed view-ordering scheme for accelerated spiral FVE (a) and its corresponding point spread function (b). In (a), each color represents a different spiral interleaf, and darker tones indicate "views" that are discarded in the partial Fourier experiments. Views aligned in *kv* are acquired sequentially throughout the cardiac cycle. Views aligned in time (same cardiac phase) are acquired in different heartbeats. In (b), each square shows the point spread function in *x*-*y* for a particular *v*-*f* coordinate. Aliasing replicas are separated from the main lobe by half of the velocity FOV, and half of the temporal frequency bandwidth, which reduces overlaps and facilitates filtering.

ROI is *s*ROI*<sup>i</sup>* (*v*, *t*) 

resolution.

the fully sampled reference.

**5.4 Acceleration results**

FVE.

interpolation (Bartels et al., 1987) along *t*.

qualitatively compared in spatial and time-velocity domains.

used as ground truth reference for a qualitative comparison.

**5.3 Acceleration experiments**

, and smoother histograms are generated by one-dimensional cubic spline

**original proposed ground truth design design reference**

The use of variable-density sampling in spiral FVE for improving spatial resolution and reducing off-resonance artifacts was evaluated in the following experiment. Three acquisitions were performed, measuring flow at the aortic valve plane of a healthy volunteer. A different spiral design was used in each acquisition (Table 2). A reduced velocity FOV (200 cm/s) was used, in order to limit each acquisition to a single feasible breath-hold. The FOV was adjusted to either the −50 to 150 cm/s range or the −150 to 50 cm/s range, depending on the flow of interest. Six *kv* encoding steps were acquired, resulting in a velocity resolution of 33 cm/s. No temporal undersampling was performed, and the acquisition was segmented across multiple heartbeats. The temporal resolution was one TR. The results were

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 21

**interleaves** 13 6 **sampling** UD VD UD **readout** 8 ms 4 ms 4 ms **field-of-view** 25 cm 6–25 cm 25 cm **resolution** 7.2 mm 3.6 mm 3.6 mm **heartbeats**† 6 18 36 **TR**‡ 12.4 ms 8.4 ms 8.4 ms

UD=uniform-density; VD=variable-density; †scan time with 6 *kv* encodes; ‡TR for 33 cm/s velocity

Table 2. Design parameters used to evaluate the use of variable-density sampling in spiral

The proposed temporal acceleration scheme was then evaluated in a second experiment. Aortic valve flow was measured using the proposed variable-density spiral design (Table 2), and the proposed view-ordering scheme (Figure 10a). The velocity resolution and FOV were set to 33 and 1200 cm/s, respectively. The data was acquired in an 18-heartbeat breath-hold, and was reconstructed using view sharing, the proposed 2D filter, and the proposed notch filter after 2D-filtering. The reduced velocity FOV results from the previous experiment were

Partial Fourier acceleration was evaluated by discarding, before reconstruction, 12 of the 36 *kv* encodes from the temporally-undersampled data from the previous experiment, as indicated in Figure 10a. The reconstructed time-velocity distribution was qualitatively compared with

Figure 12 contains results from the experiment using different spiral designs from Table 2. The data in Figure 12a was obtained using the 8 ms readout uniform-density spiral design used in section 4. The proposed variable-density design provided higher spatial resolution and reduced off-resonance artifacts, and thus better spatial localization of flow (Figure 12b). Some aliasing artifacts were observed in spatial domain (see asterisk), but these were not

Fig. 11. Aliasing in *v*-*f* space as a result of temporal undersampling. The red, green and blue ellipses illustrate the expected footprints for aortic valve flow for normal, stenotic, and regurgitant flows, respectively. Yellow dots represent peaks in the point spread function for the proposed undersampling scheme (see Figure 10). Grey ellipses represent potential aliasing components. Replicas at ±60 Hz are exact copies of the true signal. The potential aliasing at ±20 and ±40 Hz have a small footprint, because they are composed of signal from the periphery of the spatial FOV, i.e. static tissue or slow moving flow. A 2D filter (dashed lines) is capable of avoiding aliasing while preserving all signal content.

#### **5.1.3 Acceleration via partial Fourier**

Partial Fourier along the velocity dimension has been successfully used in FVE for scan time reduction, without significant loss of velocity resolution. This approach has been previously demonstrated in studies with healthy volunteers (Carvalho & Nayak, 2007; Macgowan et al., 2005) and patients (Carvalho & Nayak, 2007; Santos et al., 2007), and in phantom experiments (DiCarlo et al., 2005). Data was acquired with full coverage of *kv* space, and 33% of the data was retrospectively discarded before reconstruction (dark-colored squares in Figure 10a). The missing data was synthesized using homodyne reconstruction (Noll, Nishimura & Macovski, 1991).

#### **5.2 Data reconstruction**

Reconstruction was performed in Matlab (The MathWorks, Inc., Natick, MA, USA). The acquired data, *S*(*kx*, *ky*, *kv*, *t*), is first re-sampled onto a Cartesian grid (Jackson et al., 1991) using a Kaiser-Bessel kernel designed for the largest FOV (25 cm). Each spiral interleaf is gridded separately, and inverse Fourier transformed to form a spatial image for its corresponding *kv*-*t* coordinate, resulting in *S*(*x*, *y*, *kv*, *t*). The data corresponding to the two central *kv* values (*kv* = 0 and *kv* = <sup>1</sup> FOV*<sup>v</sup>* ) are separately filtered using a 6-tap moving average temporal filter that effectively implements view sharing (Riederer et al., 1988). A color-flow video (Riederer et al., 1991) is obtained from the filtered data. The operator draws one or multiple ROIs over the video. Pixel values within each ROI are averaged, resulting in multiple 2D datasets: *S*ROI*<sup>i</sup>* (*kv*, *t*) = ∑ROI*<sup>i</sup> <sup>x</sup>*,*<sup>y</sup> S*(*x*, *y*, *kv*, *t*). Each of these 2D datasets is filtered using the 2D filter and the notch filter described in section 5.1.2. Saturation effects (Gao et al., 1988) are compensated by normalizing the data in each cardiac phase. The data is then zero-padded along the *kv* axis, and homodyne reconstruction (Noll, Nishimura & Macovski, 1991) is used to produce each *s*ROI*<sup>i</sup>* (*v*, *t*) distribution. The time-velocity histogram for each ROI is *s*ROI*<sup>i</sup>* (*v*, *t*) , and smoother histograms are generated by one-dimensional cubic spline interpolation (Bartels et al., 1987) along *t*.

#### **5.3 Acceleration experiments**

18 Will-be-set-by-IN-TECH

stenosis

normal flow

regurgitation

−60 60

frequency (Hz)

Fig. 11. Aliasing in *v*-*f* space as a result of temporal undersampling. The red, green and blue ellipses illustrate the expected footprints for aortic valve flow for normal, stenotic, and regurgitant flows, respectively. Yellow dots represent peaks in the point spread function for the proposed undersampling scheme (see Figure 10). Grey ellipses represent potential aliasing components. Replicas at ±60 Hz are exact copies of the true signal. The potential aliasing at ±20 and ±40 Hz have a small footprint, because they are composed of signal from the periphery of the spatial FOV, i.e. static tissue or slow moving flow. A 2D filter (dashed

Partial Fourier along the velocity dimension has been successfully used in FVE for scan time reduction, without significant loss of velocity resolution. This approach has been previously demonstrated in studies with healthy volunteers (Carvalho & Nayak, 2007; Macgowan et al., 2005) and patients (Carvalho & Nayak, 2007; Santos et al., 2007), and in phantom experiments (DiCarlo et al., 2005). Data was acquired with full coverage of *kv* space, and 33% of the data was retrospectively discarded before reconstruction (dark-colored squares in Figure 10a). The missing data was synthesized using homodyne reconstruction (Noll,

Reconstruction was performed in Matlab (The MathWorks, Inc., Natick, MA, USA). The acquired data, *S*(*kx*, *ky*, *kv*, *t*), is first re-sampled onto a Cartesian grid (Jackson et al., 1991) using a Kaiser-Bessel kernel designed for the largest FOV (25 cm). Each spiral interleaf is gridded separately, and inverse Fourier transformed to form a spatial image for its corresponding *kv*-*t* coordinate, resulting in *S*(*x*, *y*, *kv*, *t*). The data corresponding to the two

temporal filter that effectively implements view sharing (Riederer et al., 1988). A color-flow video (Riederer et al., 1991) is obtained from the filtered data. The operator draws one or multiple ROIs over the video. Pixel values within each ROI are averaged, resulting in

using the 2D filter and the notch filter described in section 5.1.2. Saturation effects (Gao et al., 1988) are compensated by normalizing the data in each cardiac phase. The data is then zero-padded along the *kv* axis, and homodyne reconstruction (Noll, Nishimura & Macovski,

FOV*<sup>v</sup>* ) are separately filtered using a 6-tap moving average

(*v*, *t*) distribution. The time-velocity histogram for each

(*kv*, *t*) = ∑ROI*<sup>i</sup> <sup>x</sup>*,*<sup>y</sup> S*(*x*, *y*, *kv*, *t*). Each of these 2D datasets is filtered

velocity (cm/s)

**5.1.3 Acceleration via partial Fourier**

Nishimura & Macovski, 1991).

central *kv* values (*kv* = 0 and *kv* = <sup>1</sup>

1991) is used to produce each *s*ROI*<sup>i</sup>*

**5.2 Data reconstruction**

multiple 2D datasets: *S*ROI*<sup>i</sup>*

0

−600

<sup>0</sup> −600

lines) is capable of avoiding aliasing while preserving all signal content.

The use of variable-density sampling in spiral FVE for improving spatial resolution and reducing off-resonance artifacts was evaluated in the following experiment. Three acquisitions were performed, measuring flow at the aortic valve plane of a healthy volunteer. A different spiral design was used in each acquisition (Table 2). A reduced velocity FOV (200 cm/s) was used, in order to limit each acquisition to a single feasible breath-hold. The FOV was adjusted to either the −50 to 150 cm/s range or the −150 to 50 cm/s range, depending on the flow of interest. Six *kv* encoding steps were acquired, resulting in a velocity resolution of 33 cm/s. No temporal undersampling was performed, and the acquisition was segmented across multiple heartbeats. The temporal resolution was one TR. The results were qualitatively compared in spatial and time-velocity domains.


UD=uniform-density; VD=variable-density; †scan time with 6 *kv* encodes; ‡TR for 33 cm/s velocity resolution.

Table 2. Design parameters used to evaluate the use of variable-density sampling in spiral FVE.

The proposed temporal acceleration scheme was then evaluated in a second experiment. Aortic valve flow was measured using the proposed variable-density spiral design (Table 2), and the proposed view-ordering scheme (Figure 10a). The velocity resolution and FOV were set to 33 and 1200 cm/s, respectively. The data was acquired in an 18-heartbeat breath-hold, and was reconstructed using view sharing, the proposed 2D filter, and the proposed notch filter after 2D-filtering. The reduced velocity FOV results from the previous experiment were used as ground truth reference for a qualitative comparison.

Partial Fourier acceleration was evaluated by discarding, before reconstruction, 12 of the 36 *kv* encodes from the temporally-undersampled data from the previous experiment, as indicated in Figure 10a. The reconstructed time-velocity distribution was qualitatively compared with the fully sampled reference.

#### **5.4 Acceleration results**

Figure 12 contains results from the experiment using different spiral designs from Table 2. The data in Figure 12a was obtained using the 8 ms readout uniform-density spiral design used in section 4. The proposed variable-density design provided higher spatial resolution and reduced off-resonance artifacts, and thus better spatial localization of flow (Figure 12b). Some aliasing artifacts were observed in spatial domain (see asterisk), but these were not


Fig. 13. Temporal acceleration compared with view sharing in (left) *v*-*f* space and (right) *v*-*t* space: (a) undersampled data; (b) with 2D filtering; (c) with 2D and notch filtering (proposed approach); and (d) with view sharing (conventional approach). The 2D filter (dashed lines) removes a majority of the aliasing, and the notch filter (dotted line) removes the remaining aliasing signal (solid arrows). The proposed method removes aliasing components without noticeable loss of temporal resolution. View sharing reduces the temporal frequency bandwidth (dashed arrows), which causes temporal blurring (circles). Compare the *v*-*f*

No significant artifacts were observed when comparing the reference dataset with the 18-fold

Figure 15 presents time-velocity distributions from multiple ROIs. These distributions were reconstructed from the 18-fold undersampled dataset used in Figure 14d. Very few artifacts were observed in time-velocity histograms measured in voxels from different locations in the heart. Artifacts could be further reduced by designing different 2D filters for each ROI, based on typical characteristics of the targeted flow. For example, the artifacts observed in the descending aorta (see arrow), could be reduced by more aggressively filtering high positive-velocity components, as no ascending flow is expected in that vessel. These results, when compared with those in Figure 8, also illustrate the different improvements achieved with this approach. The spatial resolution was improved from 7 mm to 3.6 mm, and off-resonance effects were reduced. The velocity FOV was increased from ±400 to ±600 cm/s, without loss of velocity resolution (33 cm/s). The effective temporal resolution was improved from 26 ms to 9 ms, and ghosting artifacts due to view sharing were eliminated. Both

velocity (cm/s)

600 400 200 0 -200 -400 -600

−40 dB −30 dB −20 dB −10 dB 0 dB

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 23

**a**

**b**

**c**

**d**

acquisitions were performed in 12-heartbeat breath-holds.

representation in (a) with Figure 11.

accelerated result.

observed in the time-velocity distributions. A fully sampled reference is shown in Figure 12c, for comparison.

Fig. 12. Effect of variable-density sampling on image quality and spatial localization of flow: (a) original design; (b) new design; (c) ground truth reference. Top row: spatial images from the first cardiac phase; center row: time-velocity distributions measured at the aortic valve; bottom row: time-velocity distributions measured in the descending aorta. The use of higher spatial resolution and shorter readout duration improves the spatial localization of flow, which is identified by the reduced signal from static material in the time-velocity histograms (see arrows). Some aliasing artifacts were observed in spatial domain (see asterisk), but these were not observed in the time-velocity distributions.

Figure 13 contains results from the temporal acceleration experiment. Figure 13a shows the undersampled data in both *v*-*f* and *v*-*t* domains (compare this with Figure 11). Aliasing components were significantly reduced using the proposed 2D filter (dashed lines), while all of the signal energy was preserved (Figure 13b). The notch filter (dotted line) removed the majority of the remaining aliasing at ±20 and ±40 Hz (solid arrows) (Figure 13c). These results show that the proposed temporal acceleration scheme is capable of achieving 6-fold acceleration in multi-interleaf spiral FVE, without noticeable loss of temporal resolution, and without introducing significant artifacts. A result using view sharing is shown in Figure 13d, for comparison. This approach is equivalent to a moving-average low-pass filter, which reduces the temporal frequency bandwidth (dashed arrows), and causes loss of temporal resolution, perceived as blurring along *t* (circled).

Figure 14 shows a comparison between the accelerated results and the fully sampled reference. Two-fold acceleration was achieved using variable-density sampling (Figure 14b), with no noticeable artifacts in the time-velocity histogram when compared with the fully sampled reference (Figure 14a). Additional 6-fold acceleration was achieved using the proposed temporal acceleration scheme (Figure 14c). Those results were achieved in a single 18-heartbeat acquisition, while a fully sampled acquisition with the same scan parameters would require 216 heartbeats. Partial Fourier was then used to reduce the acquisition time to 12 heartbeats (i.e., by 1.5-fold), which represents a combined 18-fold acceleration (Figure 14d). 20 Will-be-set-by-IN-TECH

observed in the time-velocity distributions. A fully sampled reference is shown in Figure 12c,

500

Fig. 12. Effect of variable-density sampling on image quality and spatial localization of flow: (a) original design; (b) new design; (c) ground truth reference. Top row: spatial images from the first cardiac phase; center row: time-velocity distributions measured at the aortic valve; bottom row: time-velocity distributions measured in the descending aorta. The use of higher spatial resolution and shorter readout duration improves the spatial localization of flow, which is identified by the reduced signal from static material in the time-velocity histograms (see arrows). Some aliasing artifacts were observed in spatial domain (see asterisk), but these

Figure 13 contains results from the temporal acceleration experiment. Figure 13a shows the undersampled data in both *v*-*f* and *v*-*t* domains (compare this with Figure 11). Aliasing components were significantly reduced using the proposed 2D filter (dashed lines), while all of the signal energy was preserved (Figure 13b). The notch filter (dotted line) removed the majority of the remaining aliasing at ±20 and ±40 Hz (solid arrows) (Figure 13c). These results show that the proposed temporal acceleration scheme is capable of achieving 6-fold acceleration in multi-interleaf spiral FVE, without noticeable loss of temporal resolution, and without introducing significant artifacts. A result using view sharing is shown in Figure 13d, for comparison. This approach is equivalent to a moving-average low-pass filter, which reduces the temporal frequency bandwidth (dashed arrows), and causes loss of temporal

Figure 14 shows a comparison between the accelerated results and the fully sampled reference. Two-fold acceleration was achieved using variable-density sampling (Figure 14b), with no noticeable artifacts in the time-velocity histogram when compared with the fully sampled reference (Figure 14a). Additional 6-fold acceleration was achieved using the proposed temporal acceleration scheme (Figure 14c). Those results were achieved in a single 18-heartbeat acquisition, while a fully sampled acquisition with the same scan parameters would require 216 heartbeats. Partial Fourier was then used to reduce the acquisition time to 12 heartbeats (i.e., by 1.5-fold), which represents a combined 18-fold acceleration (Figure 14d).

time(ms) 250

**a b c**

for comparison.

velocity (cm/s)

0

−50

0

were not observed in the time-velocity distributions.

resolution, perceived as blurring along *t* (circled).

−150 0

50

150

Fig. 13. Temporal acceleration compared with view sharing in (left) *v*-*f* space and (right) *v*-*t* space: (a) undersampled data; (b) with 2D filtering; (c) with 2D and notch filtering (proposed approach); and (d) with view sharing (conventional approach). The 2D filter (dashed lines) removes a majority of the aliasing, and the notch filter (dotted line) removes the remaining aliasing signal (solid arrows). The proposed method removes aliasing components without noticeable loss of temporal resolution. View sharing reduces the temporal frequency bandwidth (dashed arrows), which causes temporal blurring (circles). Compare the *v*-*f* representation in (a) with Figure 11.

No significant artifacts were observed when comparing the reference dataset with the 18-fold accelerated result.

Figure 15 presents time-velocity distributions from multiple ROIs. These distributions were reconstructed from the 18-fold undersampled dataset used in Figure 14d. Very few artifacts were observed in time-velocity histograms measured in voxels from different locations in the heart. Artifacts could be further reduced by designing different 2D filters for each ROI, based on typical characteristics of the targeted flow. For example, the artifacts observed in the descending aorta (see arrow), could be reduced by more aggressively filtering high positive-velocity components, as no ascending flow is expected in that vessel. These results, when compared with those in Figure 8, also illustrate the different improvements achieved with this approach. The spatial resolution was improved from 7 mm to 3.6 mm, and off-resonance effects were reduced. The velocity FOV was increased from ±400 to ±600 cm/s, without loss of velocity resolution (33 cm/s). The effective temporal resolution was improved from 26 ms to 9 ms, and ghosting artifacts due to view sharing were eliminated. Both acquisitions were performed in 12-heartbeat breath-holds.

Signal-to-noise ratio was not a limiting issue for the presented application. This is in part due to the large voxel sizes and the multi-dimensional characteristics of spiral FVE. The proposed acceleration scheme reduced scan time and improved spatiotemporal resolution, and did not compromise the quality of the time-velocity histograms. Aliasing artifacts in spatial domain due to variable-density sampling are negligible in time-velocity domain (Figure 12b). The proposed *k*-*t* filters remove a majority of the temporal aliasing artifacts, and also filter high-frequency noise for high velocities (Figure 13c). Partial Fourier acceleration may cause some artifacts due to the use of a low-resolution phase estimate, thus a low acceleration factor was used (Figure 14d). Further acceleration could be achieved using parallel imaging techniques (Pruessmann et al., 2001; Samsonov et al., 2006), and further reduction in off-resonance effects could be achieved by imaging at lower field strengths.

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 25

In this chapter, we have addressed the issue of non-invasive aortic valve flow quantitation through magnetic resonance imaging. We addressed both imaging and reconstruction aspects, including accelerated acquisitions and reconstruction from undersampled data. We introduced spiral FVE, a new method for MR flow quantitation, which is capable of accurately capturing peak velocities in flow jets due to stenosis or regurgitation. Spiral FVE compared well against Doppler ultrasound, the current gold standard for cardiovascular flow imaging, and against high-resolution 2DFT phase contrast, the current MRI gold standard. Our method

Using a combination of three different techniques (variable-density spirals, temporal acceleration, and partial Fourier reconstruction), we are able to improve the spiral FVE method by 18-fold. Improvements consisted of increased velocity field-of-view, higher spatial resolution, reduced off-resonance effects, and higher temporal resolution. The improved acquisition was performed in only 12 heartbeats, whereas 216 heartbeats would be necessary to achieve such improvements without acceleration. No significant artifacts were observed. Magnetic resonance imaging is potentially the most appropriate technique for addressing all aspects of cardiovascular disease examination. The evaluation of valvular disease and intracardiac flow will be a necessary capability in a comprehensive cardiac MR examination. The imaging and reconstruction techniques proposed in this chapter can be an important

Ahn, C. B., Park, J. H. & Cho, Z. H. (1986). High-speed spiral-scan echo planar NMR imaging,

Bartels, R. H., Beatty, J. C. & Barsky, B. A. (1987). *An introduction to splines for use in computer graphics and geometric modelling*, Morgan Kaufmann, San Francisco, California. Bernstein, M. A., Fain, S. B. & Riederer, S. J. (2001). Effect of windowing and zero-filled

Bernstein, M. A., King, K. F. & Zhou, X. J. (2004). *Handbook of MRI Pulse Sequences*, Academic

Carvalho, J. L. A., DiCarlo, J. C., Kerr, A. B. & Nayak, K. S. (2007). Reconstruction

reconstruction of MRI data on spatial resolution and acquisition strategy, *J Magn*

of variable-density data in Fourier velocity encoding, *Proc, ISMRM, 15th Annual*

was demonstrated in both healthy volunteers and in a patient.

contribution towards making such exam feasible.

*IEEE Trans Med Imaging* 5(1): 2–7.

*Reson Imaging* 14(3): 270–280.

*Meeting*, Berlin, p. 2514.

**6. Conclusion**

**7. References**

Press.

Fig. 14. Accelerated spiral FVE results: (a) fully sampled reference (36 heartbeats); (b) 2-fold acceleration, using variable-density sampling (18 heartbeats); (c) 12-fold acceleration, using variable-density sampling and temporal acceleration (18 heartbeats); and (d) 18-fold acceleration, using variable-density sampling, temporal acceleration, and partial Fourier (12 heartbeats). A reduced velocity FOV (200 cm/s) was used in (a) and (b) to limit the acquisition to a single feasible breath-hold. If acquiring a full 1200 cm/s velocity FOV, as in (c) and (d), the total acquisition time for (a) and (b) would have been 216 and 108 heartbeats, respectively. All other scan parameters were identical for the four acquisitions. No significant differences were observed when comparing the 18-fold accelerated result with the fully sampled reference.

Fig. 15. Flow in multiple ROIs measured in a single 12-heartbeat spiral FVE acquisition using 18-fold acceleration. Voxels of interest are highlighted. Time-velocity histograms reveal only minimal artifacts (see arrow). The data illustrates significant improvements in spatial resolution, velocity field-of-view, and temporal resolution (compare with Figure 8).

Signal-to-noise ratio was not a limiting issue for the presented application. This is in part due to the large voxel sizes and the multi-dimensional characteristics of spiral FVE. The proposed acceleration scheme reduced scan time and improved spatiotemporal resolution, and did not compromise the quality of the time-velocity histograms. Aliasing artifacts in spatial domain due to variable-density sampling are negligible in time-velocity domain (Figure 12b). The proposed *k*-*t* filters remove a majority of the temporal aliasing artifacts, and also filter high-frequency noise for high velocities (Figure 13c). Partial Fourier acceleration may cause some artifacts due to the use of a low-resolution phase estimate, thus a low acceleration factor was used (Figure 14d). Further acceleration could be achieved using parallel imaging techniques (Pruessmann et al., 2001; Samsonov et al., 2006), and further reduction in off-resonance effects could be achieved by imaging at lower field strengths.

## **6. Conclusion**

22 Will-be-set-by-IN-TECH

Fig. 14. Accelerated spiral FVE results: (a) fully sampled reference (36 heartbeats); (b) 2-fold acceleration, using variable-density sampling (18 heartbeats); (c) 12-fold acceleration, using variable-density sampling and temporal acceleration (18 heartbeats); and (d) 18-fold

acceleration, using variable-density sampling, temporal acceleration, and partial Fourier (12

acquisition to a single feasible breath-hold. If acquiring a full 1200 cm/s velocity FOV, as in (c) and (d), the total acquisition time for (a) and (b) would have been 216 and 108 heartbeats,

significant differences were observed when comparing the 18-fold accelerated result with the

Fig. 15. Flow in multiple ROIs measured in a single 12-heartbeat spiral FVE acquisition using 18-fold acceleration. Voxels of interest are highlighted. Time-velocity histograms reveal only minimal artifacts (see arrow). The data illustrates significant improvements in spatial resolution, velocity field-of-view, and temporal resolution (compare with Figure 8).

velocity (cm/s)

0

600

velocity (cm/s)

0

600

time (ms) 0 400 800 −600

0 400 800 −600

heartbeats). A reduced velocity FOV (200 cm/s) was used in (a) and (b) to limit the

respectively. All other scan parameters were identical for the four acquisitions. No

**abcd**

**0 400 800 time (ms)**

**600**

**velocity (cm/s)**

**-600**

fully sampled reference.

velocity (cm/s)

velocity (cm/s)

0

600

0

600

0 400 800 −600

time (ms) 0 400 800 −600

In this chapter, we have addressed the issue of non-invasive aortic valve flow quantitation through magnetic resonance imaging. We addressed both imaging and reconstruction aspects, including accelerated acquisitions and reconstruction from undersampled data. We introduced spiral FVE, a new method for MR flow quantitation, which is capable of accurately capturing peak velocities in flow jets due to stenosis or regurgitation. Spiral FVE compared well against Doppler ultrasound, the current gold standard for cardiovascular flow imaging, and against high-resolution 2DFT phase contrast, the current MRI gold standard. Our method was demonstrated in both healthy volunteers and in a patient.

Using a combination of three different techniques (variable-density spirals, temporal acceleration, and partial Fourier reconstruction), we are able to improve the spiral FVE method by 18-fold. Improvements consisted of increased velocity field-of-view, higher spatial resolution, reduced off-resonance effects, and higher temporal resolution. The improved acquisition was performed in only 12 heartbeats, whereas 216 heartbeats would be necessary to achieve such improvements without acceleration. No significant artifacts were observed.

Magnetic resonance imaging is potentially the most appropriate technique for addressing all aspects of cardiovascular disease examination. The evaluation of valvular disease and intracardiac flow will be a necessary capability in a comprehensive cardiac MR examination. The imaging and reconstruction techniques proposed in this chapter can be an important contribution towards making such exam feasible.

#### **7. References**


Irarrazabal, P., Hu, B. S., Pauly, J. M. & Nishimura, D. G. (1993). Spatially resolved and localized real-time velocity distribution, *Magn Reson Med* 30(2): 207–212. Jackson, J. I., Meyer, C. H., Nishimura, D. G. & Macovski, A. (1991). Selection of a convolution

Rapid Quantitation of Aortic Valve Flow Using Spiral Fourier Velocity Encoded MRI 27

Liu, C.-Y., Varadarajan, P., Pohost, G. M. & Nayak, K. S. (2008). Real-time color-flow MRI at

Macgowan, C. K., Kellenberger, C. J., Detsky, J. S., Roman, K. & Yoo, S. J. (2005). Real-time

Madore, B., Glover, G. H. & Pelc, N. J. (1999). Unaliasing by Fourier-encoding the overlaps

Meyer, C. H., Hu, B. S., Nishimura, D. G. & Macovski, A. (1992). Fast spiral coronary artery

Mohiaddin, R. H., Gatehouse, P. D., Henien, M. & Firmin, D. N. (1997). Cine MR

Moran, P. R. (1982). A flow velocity zeugmatographic interlace for NMR imaging in humans,

Moran, P. R., Moran, R. A. & Karstaedt, N. (1985). Verification and evaluation of internal

Nayak, K. S., Pauly, J. M., Kerr, A. B., Hu, B. S. & Nishimura, D. G. (2000). Real-time color flow

Nayler, G. L., Firmin, D. N. & Longmore, D. B. (1986). Blood flow imaging by cine magnetic

Nishimura, D. G. (2010). *Principles of Magnetic Resonance Imaging*, Stanford University, Palo

Noll, D. C., Meyer, C. H., Pauly, J. M., Nishimura, D. G. & Macovski, A. (1991). A homogeneity

Noll, D. C., Nishimura, D. G. & Macovski, A. (1991). Homodyne detection in magnetic

O'Donnell, M. (1985). NMR blood flow imaging using multiecho, phase contrast sequences,

Pruessmann, K. P., Weiger, M., Bornert, P. & Boesiger, P. (2001). Advances in sensitivity encoding with arbitrary k-space trajectories, *Magn Reson Med* 46(4): 638–651. Rebergen, S. A., van der Wall, E. E., Doornbos, J. & de Roos, A. (1993). Magnetic resonance

Riederer, S. J., Tasciyan, T., Farzaneh, F., Lee, J. N., Wright, R. C. & Herfkens, R. J. (1988). MR

Riederer, S. J., Wright, R. C., Ehman, R. L., Rossman, P. J., Holsinger-Bampton, A. E.,

Samsonov, A. A., Block, W. F., Arunachalam, A. & Field, A. S. (2006). Advances in locally constrained k-space-based parallel MRI, *Magn Reson Med* 55(2): 431–438.

resonance imaging, *IEEE Trans Med Imaging* 10(2): 154–163.

fluoroscopy: technical feasibility, *Magn Reson Med* 8(1): 1–15.

correction method for magnetic resonance imaging with time-varying gradients,

measurement of velocity and flow: technique, validation, and cardiovascular

Hangiangdreu, N. J. & Grimm, R. C. (1991). Real-time interactive color flow MR

*Magn Reson Med* 42(5): 813–828.

*Magn Reson Imaging* 1(4): 197–203.

method, *Radiology* 154(2): 433–441.

MRI, *Magn Reson Med* 43: 251–258.

*IEEE Trans Med Imaging* 10(4): 629–637.

applications, *Am Heart J* 126(6): 1439–1456.

imaging, *Radiology* 181: 33–39.

Alto, California.

*Med Phys* 12(1): 59–64.

resonance, *J Comput Assist Tomogr* 10: 715–722.

imaging, *Magn Reson Med* 28(2): 202–213.

echocardiography., *J Magn Reson Imaging* 7(4): 657–663.

function for Fourier inversion using gridding, *IEEE Trans Med Imaging* 10(3): 473–478.

3 T using variable-density spiral phase contrast, *Magn Reson Imaging* 26(5): 661–666.

Fourier velocity encoding: an in vivo evaluation, *J Magn Reson Imaging* 21(3): 297–304.

using the temporal dimension (UNFOLD), applied to cardiac imaging and fMRI,

Fourier velocimetry of blood flow through cardiac valves: comparison with Doppler

flow and motion. true magnetic resonance imaging by the phase gradient modulation


24 Will-be-set-by-IN-TECH

Carvalho, J. L. A. & Nayak, K. S. (2007). Rapid quantitation of cardiovascular flow using

Carvalho, J. L. A., Nielsen, J. F. & Nayak, K. S. (2010). Feasibility of in vivo measurement of

Chia, J. M., Fischer, S. E., Wickline, S. A. & Lorenz, C. H. (2000). Arrhythmia rejection using a VCG-based triggering algorithm, *Proc, ISMRM, 8th Annual Meeting*, Denver, p. 201. Clarke, G. D., Hundley, W. G., McColl, R. W., Eckels, R., Smith, D., Chaney, C., Li, H. F. &

DiCarlo, J. C., Hargreaves, B. A., Nayak, K. S., Hu, B. S., Pauly, J. M. & Nishimura, D. G. (2005).

Dumoulin, C. L., Souza, S. P., Hardy, C. J. & Ash, S. A. (1991). Quantitative measurement of

Feinberg, D. A., Crooks, L. E., Sheldon, P., Hoenninger III, J., Watts, J. & Arakawa, M. (1985).

Feinberg, D. A. & Oshio, K. (1992). Gradient-echo shifting in fast MRI techniques (GRASE

Galea, D., Lauzon, M. L. & Drangova, M. (2002). Peak velocity determination using fast

Gamper, U., Boesiger, P. & Kozerke, S. (2008). Compressed sensing in dynamic MRI, *Magn*

Gao, J. H., Holland, S. K. & Gore, J. C. (1988). Nuclear magnetic resonance signal from flowing nuclei in rapid imaging using gradient echoes, *Med Phys* 15(6): 809–814. Glover, G. H. & Pelc, N. J. (1988). A rapid-gated cine MRI technique, *Magn Reson Annu*

Hahn, E. L. (1960). Detection of sea-water motion by nuclear precession, *J Geophys Res*

Hansen, M. S., Baltes, C., Tsao, J., Kozerke, S., Pruessmann, K. P., Boesiger, P. & Pedersen,

Hargreaves, B. A. (2001). *Spin-Manipulation Methods for Efficient Magnetic Resonance Imaging*,

Hennig, J., Mueri, M., Brunner, P. & Friedburg, H. (1988). Fast and exact flow measurements

Holsinger, A. E., Wright, R. C., Riederer, S. J., Farzaneh, F., Grimm, R. C. & Maier, J. K. (1990). Real-time interactive magnetic resonance imaging, *Magn Reson Med* 14(3): 547–553. Hoskins, P. R. (1996). Accuracy of maximum velocity estimates made using Doppler

Hu, B. S., Pauly, J. M. & Macovski, A. (1993). Localized real-time velocity spectra

with the fast Fourier flow technique, *Magn Reson Med* 6(4): 369–372.

velocity-spatio-temporal correlations, *MAGMA* 17(2): 86–94.

ultrasound systems, *Br J Radiol* 69(818): 172–177.

determination, *Magn Reson Med* 30(3): 393–398.

E. M. (2004). Accelerated dynamic Fourier velocity encoding by exploiting

57(4): 639–646.

63(6): 1537–1547.

21(2): 242–250.

*Reson Med* 2(6): 555–566.

*Reson Med* 59(2): 365–373.

PhD thesis, Stanford University.

*Reson* 97: 177–183.

pp. 299–333.

65(2): 776–777.

*J Magn Reson Imaging* 6(5): 733–742.

slice-selective Fourier velocity encoding with spiral readouts, *Magn Reson Med*

carotid wall shear rate using spiral Fourier velocity encoded MRI, *Magn Reson Med*

Peshock, R. M. (1996). Velocity-encoded, phase-difference cine MRI measurements of coronary artery flow: dependence of flow accuracy on the number of cine frames,

Variable-density one-shot Fourier velocity encoding, *Magn Reson Med* 54(3): 645–655.

blood flow using cylindrically localized Fourier velocity encoding, *Magn Reson Med*

Magnetic resonance imaging the velocity vector components of fluid flow, *Magn*

imaging) for correction of field inhomogeneity errors and chemical shift, *J Magn*

Fourier velocity encoding with minimal spatial encoding, *Med Phys* 29(8): 1719–1728.


**2** 

*Belgium* 

**State-Of-The-Art Methods for the** 

*2Department of Mechanics, University College Ghent,* 

*3IBiTech-bioMMeda, Ghent University,* 

**Numerical Simulation of Aortic BMHVs** 

Sebastiaan Annerel1, Tom Claessens2, Peter Van Ransbeeck2, Patrick Segers3, Pascal Verdonck3 and Jan Vierendeels1

*1Department of Flow, Heat and Combustion Mechanics, Ghent University,* 

Since the first clinical implantation of an artificial aortic valve by Dr. Charles A. Hufnagel in 1952 (Hufnagel et al., 1954), the use of such prostheses has gained strong interest and has become a routine treatment for severe heart valve failure. During the past 60 years, various mechanical heart valve designs have been developed for use in both aortic and mitral positions (Butany et al., 2003; Aslam et al., 2007). Nowadays, bileaflet mechanical heart valves (BMHVs) are widely preferred for aortic valve replacement because of their long lifespan. However, current BMHVs still induce pannus and thromboembolism, among other undesired side effects, which are believed to be due to non-physiological flow and

One way to gain insight into the dynamics of a BMHV in order to improve its design is by experimental testing (Grigioni et al., 2004). Usually, *in vitro* testing is used, in which the functioning of the valve is assessed, for example, by using Doppler echocardiography (Dumont et al., 2002; Verdonck et al., 2002) or by visualizing the temporal and spatial flow field through velocimetry, like the laser Doppler anemometry (LDA) technique (Browne et al., 2000; Akutsu et al., 2001) or the particle image velocimetry (PIV) technique (Browne et al., 2000; Kaminsky et al., 2007). Also, the spectrum of the valve noise can be analyzed, as is done, for example, in Masson & Rieu (1998). Experimental *in vivo* testing is another option, using echocardiography and Doppler ultrasound to investigate the behavior of the valve after implantation in human patients (Bech-Hanssen, 2001; Aslam et al., 2007; Aljassim et al.,

Numerical ("*in silico"*) methods can provide an alternative way to obtain relevant and detailed information for valve design optimization, since they are capable of solving the valve dynamics with a high degree of resolution in time and space (Kelly et al., 1999; Grigioni et al., 2004; Yoganathan et al., 2005; Dasi et al., 2009; Sotiropoulous & Borazjani, 2009). Moreover, they are considerably less time-consuming and less expensive during the research and development phase compared with experimental testing (Dasi et al., 2009) and are, therefore, particularly efficient for sensitivity studies (Verdonck, 2002). Unfortunately, the numerical simulation of a BMHV is a complex fluid-structure interaction (FSI) problem

turbulence generated by the valve leaflets (Sotiropoulous & Borazjani, 2009).

2008; Zogbi et al., 2009) or in animals (Yin et al., 2006).

**1. Introduction** 


## **State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs**

Sebastiaan Annerel1, Tom Claessens2, Peter Van Ransbeeck2, Patrick Segers3, Pascal Verdonck3 and Jan Vierendeels1 *1Department of Flow, Heat and Combustion Mechanics, Ghent University, 2Department of Mechanics, University College Ghent, 3IBiTech-bioMMeda, Ghent University, Belgium* 

## **1. Introduction**

26 Will-be-set-by-IN-TECH

28 Aortic Valve

Santos, J. M., Kerr, A. B., Lee, D., McConnell, M. V., Yang, P. C., Hu, B. S. & Pauly, J. M. (2007).

Singer, J. R. (1959). Blood flow rates by nuclear magnetic resonance measurements, *Science*

Singer, J. R. & Crooks, L. E. (1983). Nuclear magnetic resonance blood flow measurements in

Tang, C., Blatter, D. D. & Parker, D. L. (1993). Accuracy of phase-contrast flow measurements in the presence of partial-volume effects, *J Magn Reson Imaging* 3(3): 377–385. Tsai, C.-M. & Nishimura, D. G. (2000). Reduced aliasing artifacts using variable-density

Tsai, C.-M., Olcott, E. W. & Nishimura, D. G. (1999). Flow quantification using

van Dijk, P. (1984). Direct cardiac NMR imaging of heart wall and blood flow velocity, *J Comput*

Winkler, A. J. & Wu, J. (1995). Correction of intrinsic spectral broadening errors in Doppler

Zur, Y., Wood, M. L. & Neuringer, L. J. (1991). Spoiling of transverse magnetization in

low-spatial-resolution and low-velocity-resolution velocity images, *Magn Reson Med*

peak velocity measurements made with phased sector and linear array transducers,

p. 2551.

130(3389): 1652–1653.

42(4): 682–690.

*Assist Tomogr* 8(3): 429–436.

*Ultrasound Med Biol* 21(8): 1029–1035.

the human brain, *Science* 221(4611): 654–656.

k-space sampling trajectories, *Magn Reson Med* 43: 452–458.

Tsao, J. (2002). On the UNFOLD method, *Magn Reson Med* 47(1): 202–207.

steady-state sequences, *Magn Reson Med* 21(2): 251–263.

Comprehensive valve evaluation system, *Proc, ISMRM, 15th Annual Meeting*, Berlin,

Since the first clinical implantation of an artificial aortic valve by Dr. Charles A. Hufnagel in 1952 (Hufnagel et al., 1954), the use of such prostheses has gained strong interest and has become a routine treatment for severe heart valve failure. During the past 60 years, various mechanical heart valve designs have been developed for use in both aortic and mitral positions (Butany et al., 2003; Aslam et al., 2007). Nowadays, bileaflet mechanical heart valves (BMHVs) are widely preferred for aortic valve replacement because of their long lifespan. However, current BMHVs still induce pannus and thromboembolism, among other undesired side effects, which are believed to be due to non-physiological flow and turbulence generated by the valve leaflets (Sotiropoulous & Borazjani, 2009).

One way to gain insight into the dynamics of a BMHV in order to improve its design is by experimental testing (Grigioni et al., 2004). Usually, *in vitro* testing is used, in which the functioning of the valve is assessed, for example, by using Doppler echocardiography (Dumont et al., 2002; Verdonck et al., 2002) or by visualizing the temporal and spatial flow field through velocimetry, like the laser Doppler anemometry (LDA) technique (Browne et al., 2000; Akutsu et al., 2001) or the particle image velocimetry (PIV) technique (Browne et al., 2000; Kaminsky et al., 2007). Also, the spectrum of the valve noise can be analyzed, as is done, for example, in Masson & Rieu (1998). Experimental *in vivo* testing is another option, using echocardiography and Doppler ultrasound to investigate the behavior of the valve after implantation in human patients (Bech-Hanssen, 2001; Aslam et al., 2007; Aljassim et al., 2008; Zogbi et al., 2009) or in animals (Yin et al., 2006).

Numerical ("*in silico"*) methods can provide an alternative way to obtain relevant and detailed information for valve design optimization, since they are capable of solving the valve dynamics with a high degree of resolution in time and space (Kelly et al., 1999; Grigioni et al., 2004; Yoganathan et al., 2005; Dasi et al., 2009; Sotiropoulous & Borazjani, 2009). Moreover, they are considerably less time-consuming and less expensive during the research and development phase compared with experimental testing (Dasi et al., 2009) and are, therefore, particularly efficient for sensitivity studies (Verdonck, 2002). Unfortunately, the numerical simulation of a BMHV is a complex fluid-structure interaction (FSI) problem

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 31

Secondly, the blood flow through the valve is calculated, which is governed by the conservation of mass and the Navier-Stokes equations. For the unsteady flow of an

( ) *<sup>v</sup>* <sup>2</sup> *vv p v g <sup>t</sup>*

G

time, *μ* = dynamic viscosity and *p* = pressure. In the case of blood flow through a BMHV, it is usually assumed that the hinge axis is in the direction of gravity. Thus, the gravity has no effect on the moment about the hinge axis (Sotiropoulos et al., 2009) and, therefore, its

Computational fluid dynamics (CFD) are used to solve Equation (2) for the entire fluid domain. From the resulting flow field, the pressure and viscous forces at the fluid-structure interface are derived. These forces are integrated at the fluid-structure interface giving the

Finally, the (kinematic and dynamic) equilibrium equations at the fluid-structure interface need to be solved. The kinematic equilibrium states that the angular position of the fluid at the fluid-structure interface (i.e. *θf,i*) should be equal to the angular position of the structural

> ,1 ,1 ,2 ,2 *f s f s*

Dynamic equilibrium also needs to be achieved. When the hinges of the valve are modeled as frictionless, the structural moment *Ms,i* acting on each leaflet *i* should be equal to the

θ

θ

,1 ,1 ,2 ,2 *f s f s*

,1 1 ,1 ,2 2 ,2 *f s f s*

θ

θ

*M I M I*

⎨

<sup>⎧</sup> = ⋅ <sup>⎪</sup>

<sup>⎪</sup> = ⋅ <sup>⎩</sup>

Both equations of equilibrium at the fluid-structure interface are solved using FSI methods. The subscripts *s* and *f* in Equation (5) are left out from here on. In the following, the pressure

> 1 11 2 22

θ

θ

*M I M I*

⎪⎧ <sup>=</sup> <sup>⋅</sup> <sup>⎨</sup> ⎪ = ⋅ ⎩

*M M M M* ⎧ = ⎪ <sup>⎨</sup> <sup>=</sup> ⎪⎩

θ

⎧ = ⎪ <sup>⎨</sup> <sup>=</sup> ⎪⎩

θ

pressure and viscous moment exerted by the flow, indicated by *Mf,i*:

and viscous moment *Mf,i* and the structural acceleration

referred to as "the moment *Mi*" and "the angular acceleration *<sup>i</sup>*

<sup>∂</sup> + ∇ ⋅ = −∇ + ∇ +

 μρ

<sup>G</sup> (2a)

= flow velocity vector, *ρ* = fluid density, *t* =

(3)

(4)

(5)

(6)

*<sup>s</sup>*,1 will thus respectively be

". With this change in

θ

θ

<sup>G</sup> GG GG (2b)

incompressible fluid, the differential equations to be solved are given by:

ρ

pressure and viscous moment *Mf,i* about the hinge axis acting on the interface.

0 ∇ ⋅ = *v*

(Hirt et al., 1997; Donea et al., 2004), in which *v*

influence (last term in Equation (2b)) is neglected.

leaflet at the interface (i.e. *θs,i*):

or with Equation (1):

notation, Equation (5) becomes:

ρ

∂

because the movement of the leaflets strongly interacts with the surrounding fluid motion; therefore, the equilibrium at the fluid-structure interface needs to be taken into account. In this chapter, a review of numerical FSI methods for BMHVs is given. Subsequently, the general dynamics and flow fields of BMHVs are discussed and illustrated by numerical simulations. This flow field typically consists of three jets. Furthermore, the design optimization challenges are described. High-flow-velocity gradients give rise to high shear stresses that can induce blood damage (Yoganathan et al., 2004). Therefore, the blood damage is discussed. The flow through the hinge region is of special interest. Finally, the cavitation phenomenon in BMHVs is discussed, because it can induce blood damage as well as structural failure (due to pitting and erosion).

## **2. A review of FSI methods to simulate the dynamics of a BMHV**

When numerically simulating a BMHV, three problems need to be solved, namely the structural problem, the flow problem, and the interaction of the fluid and the structure at the fluid-structure interface. In the following, each of these three problems is discussed in detail.

Since the leaflets of a BMHV have a small moment of inertia and are very stiff, they are usually assumed to be rigid. A BMHV can thus be modeled as a rigid casing in which two separate rigid leaflets rotate around their hinge axes (see Fig. 1). Because the position of each rigid leaflet is solely determined by its opening angle, the bileaflet valve has two degrees of freedom.

Fig. 1. View on the ATS Open PivotTM Standard Heart Valve with leaflets (marked in black) in the open position. The casing is visible (in white) with the blocking mechanism at the hinges

The movement of a rigid leaflet *i* is governed by Newton's Second Law, which states that the (structural) moment *Ms,i* about its rotation axis must be in equilibrium with the product of its moment of inertia *Ii* and its angular acceleration *s i*, θ . For two leaflets, this results in the following two equations:

$$\begin{cases} M\_{s,1} = I\_1 \cdot \ddot{\theta}\_{s,1} \\ M\_{s,2} = I\_2 \cdot \ddot{\theta}\_{s,2} \end{cases} \tag{1}$$

Secondly, the blood flow through the valve is calculated, which is governed by the conservation of mass and the Navier-Stokes equations. For the unsteady flow of an incompressible fluid, the differential equations to be solved are given by:

$$\nabla \cdot \vec{\upsilon} = 0 \tag{2a}$$

$$
\rho \frac{\partial \vec{v}}{\partial t} + \rho \nabla \left( \vec{v} \vec{v} \right) = -\nabla p + \mu \nabla^2 \vec{v} + \rho \vec{g} \tag{2b}
$$

(Hirt et al., 1997; Donea et al., 2004), in which *v* G = flow velocity vector, *ρ* = fluid density, *t* = time, *μ* = dynamic viscosity and *p* = pressure. In the case of blood flow through a BMHV, it is usually assumed that the hinge axis is in the direction of gravity. Thus, the gravity has no effect on the moment about the hinge axis (Sotiropoulos et al., 2009) and, therefore, its influence (last term in Equation (2b)) is neglected.

Computational fluid dynamics (CFD) are used to solve Equation (2) for the entire fluid domain. From the resulting flow field, the pressure and viscous forces at the fluid-structure interface are derived. These forces are integrated at the fluid-structure interface giving the pressure and viscous moment *Mf,i* about the hinge axis acting on the interface.

Finally, the (kinematic and dynamic) equilibrium equations at the fluid-structure interface need to be solved. The kinematic equilibrium states that the angular position of the fluid at the fluid-structure interface (i.e. *θf,i*) should be equal to the angular position of the structural leaflet at the interface (i.e. *θs,i*):

$$\begin{cases} \theta\_{f,1} = \theta\_{s,1} \\ \theta\_{f,2} = \theta\_{s,2} \end{cases} \tag{3}$$

Dynamic equilibrium also needs to be achieved. When the hinges of the valve are modeled as frictionless, the structural moment *Ms,i* acting on each leaflet *i* should be equal to the pressure and viscous moment exerted by the flow, indicated by *Mf,i*:

$$\begin{cases} M\_{f,1} = M\_{s,1} \\ M\_{f,2} = M\_{s,2} \end{cases} \tag{4}$$

or with Equation (1):

30 Aortic Valve

because the movement of the leaflets strongly interacts with the surrounding fluid motion; therefore, the equilibrium at the fluid-structure interface needs to be taken into account. In this chapter, a review of numerical FSI methods for BMHVs is given. Subsequently, the general dynamics and flow fields of BMHVs are discussed and illustrated by numerical simulations. This flow field typically consists of three jets. Furthermore, the design optimization challenges are described. High-flow-velocity gradients give rise to high shear stresses that can induce blood damage (Yoganathan et al., 2004). Therefore, the blood damage is discussed. The flow through the hinge region is of special interest. Finally, the cavitation phenomenon in BMHVs is discussed, because it can induce blood damage as well

When numerically simulating a BMHV, three problems need to be solved, namely the structural problem, the flow problem, and the interaction of the fluid and the structure at the fluid-structure interface. In the following, each of these three problems is discussed in

Since the leaflets of a BMHV have a small moment of inertia and are very stiff, they are usually assumed to be rigid. A BMHV can thus be modeled as a rigid casing in which two separate rigid leaflets rotate around their hinge axes (see Fig. 1). Because the position of each rigid leaflet is solely determined by its opening angle, the bileaflet valve has two degrees of

Fig. 1. View on the ATS Open PivotTM Standard Heart Valve with leaflets (marked in black) in the open position. The casing is visible (in white) with the blocking mechanism at the

The movement of a rigid leaflet *i* is governed by Newton's Second Law, which states that the (structural) moment *Ms,i* about its rotation axis must be in equilibrium with the product

> ,1 1 ,1 ,2 2 ,2 *s s s s*

θ

θ

*M I M I*

⎨

⎧ = ⋅ ⎪

⎪ = ⋅ ⎩

θ

. For two leaflets, this results in the

(1)

of its moment of inertia *Ii* and its angular acceleration *s i*,

**2. A review of FSI methods to simulate the dynamics of a BMHV** 

as structural failure (due to pitting and erosion).

detail.

freedom.

hinges

following two equations:

$$\begin{cases} M\_{f,1} = I\_1 \cdot \ddot{\theta}\_{s,1} \\ M\_{f,2} = I\_2 \cdot \ddot{\theta}\_{s,2} \end{cases} \tag{5}$$

Both equations of equilibrium at the fluid-structure interface are solved using FSI methods. The subscripts *s* and *f* in Equation (5) are left out from here on. In the following, the pressure and viscous moment *Mf,i* and the structural acceleration θ*<sup>s</sup>*,1 will thus respectively be referred to as "the moment *Mi*" and "the angular acceleration *<sup>i</sup>* θ ". With this change in notation, Equation (5) becomes:

$$\begin{cases} M\_1 = I\_1 \cdot \ddot{\theta}\_1 \\ M\_2 = I\_2 \cdot \ddot{\theta}\_2 \end{cases} \tag{6}$$

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 33

the ALE approach is its accuracy, because the grid is aligned with the fluid-structure interface. However, the use of remeshing (and thus interpolation) introduces artificial diffusivity and can become expensive for complex three-dimensional geometries. Several studies have used the ALE approach to simulate the dynamics of the ATS Open PivotTM Standard Heart Valve (Dumont et al., 2005, 2007), the St. Jude MedicalTM BMHV (Penrose et al., 2002; Redaelli et al., 2004; Dumont et al., 2007; Guivier et al., 2007, 2009; Nobili et al., 2007, 2008; Choi et al., 2009; Hong et al., 2009), and other valve types (Makhijani et al., 1997;

Secondly, one can classify each FSI simulation by using a partitioned solver or a monolithic solver. In the monolithic approach, the entire FSI problem is simultaneously simulated by

This is in contrast to the partitioned approach, which solves the flow and the structural problem separately and, therefore, mostly uses different specialized solvers. The partitioned approach is used to simulate heart valves in Makhijani et al. (1997), Penrose et al. (2002), Redaelli et al. (2004), Dumont et al. (2005, 2007), Vierendeels et al. (2005, 2007), Bang et al. (2006), Guivier et al. (2007, 2009), Nobili et al. (2007, 2008), Morsi et al. (2007), Tai et al. (2007), Borazjani et al. (2008), Diniz dos Santos et al. (2008), Astorino et al. (2009), Choi et al. (2009), De Tullio et al. (2009), Hong et al. (2009) and, finally, Xia et al. (2009). In order to obtain the interaction between the fluid and the structure, data exchange at the fluid-structure interface and a coupling scheme between the separated solvers are needed. Unfortunately, not every coupling scheme converges quickly. The instability of coupling schemes without relaxation is analytically explained in Vierendeels et al. (2005) and Borazjani et al. (2008) for the case of BMHVs, and it is also demonstrated by the flow through arteries (Causin et al., 2005; Degroote et al., 2008, 2010). It is concluded that relaxation can be used to obtain stable and efficient approximations for the subsequent angular accelerations of the valve leaflets. Several coupling schemes with relaxation have thus been developed, and they can be divided into loose and strong

In the loose coupling methods, only one coupling iteration is needed in each time-step, since the solution of the flow field at time-step *n* is used to calculate the angular accelerations of

( ) 11 1 <sup>1</sup>

 ω

*n n nn i i i ii*

is unstable (Causin et al., 2005; Vierendeels et al., 2005; Borazjani et al., 2008), as mentioned above. The loose coupling formulation is often used to simulate heart valves (Redaelli et al., 2004; Morsi et al., 2007; Nobili et al., 2007; Tai et al., 2007; Borazjani et al., 2008; Xia et al., 2009). It has the main benefit of a low computational cost because of the lack of a coupling iteration loop within each time-step. However, this lack implies that dynamic equilibrium at the fluid-structure interface (Equation (6)) is not necessarily achieved, which leads to unphysical oscillations in the leaflet movement and the flow and pressure field, as described

 θω

<sup>+</sup> < 1) is necessary, since the scheme without relaxation ( *<sup>n</sup>* <sup>1</sup>

*n*

*M I*

++ + =− ⋅ + ⋅ (8)

*i*

*i* ω<sup>+</sup> = 1)

Vierendeels et al., 2005, 2007; Bang et al., 2006; Morsi et al., 2007).

**2.2 Monolithic solver versus partitioned solver** 

one solver.

coupling schemes.

A relaxation factor ( *<sup>n</sup>* <sup>1</sup>

in Annerel et al. (2011).

the leaflets for the next time-step *n*+1:

*i* ω

θ

In the remainder of this section, classifications of FSI methods used in literature to simulate a BMHV will be made, and the characteristics of each of the mentioned methods will be explained.

#### **2.1 Fixed grid techniques versus moving grid techniques**

A first classification concerns the kinematical description of the domain. For the structural problem, the motion is usually described by the Lagrangian method, where the computational grid moves with the material velocity. This is in contrast to a fluid domain, in which the motion is generally described by the Eulerian method and, therefore, the computational grid does not deform. In case of FSI, both methods can be combined in several approaches in order to describe the motion of the domain.

One approach is the "fixed grid" method, in which an immersed structure is allowed to "fictitiously" move through the Eulerian fluid grid in a Lagrangian way. The influence of the structure boundary on the fluid outside the structure is calculated by introducting body force sources into the Navier-Stokes equations, while keeping the velocity of the fictitiously overlapped fluid coupled to the structural velocity. Since the fluid grid is fixed, there is no need for remeshing and grid adaption. However, because the fluid-structure interface is not necessarily aligned with the spatial discretization and the data are, as such, interpolated, the flow field (and thus shear stresses) at this interface is not accurately calculated. Several fixed grid methods have been developed, such as the immersed boundary (IB) method, first proposed by Peskin (1972) for the simulation of heart valves. Borazjani et al. (2008) used the sharp interface CURVIB-method for simulating a BMHV. Other IB techniques were developed and used by Tai et al. (2007), De Tullio et al. (2009) and Xia et al. (2009). Also, the fictitious domain (FD) method can be used to simulate flexible heart valves (De Hart et al., 2000, 2003; Diniz dos Santos et al., 2008; Astorino et al., 2009). This fixed grid method uses Lagrange multipliers to impose the kinematic constraints. Van Loon et al. (2004) improved the accuracy of the FD method at the fluid-structure interface by proposing a local fluid grid adaption at the structure boundary.

Another approach is to use the arbitrary Lagrangian-Eulerian (ALE) method for the fluid domain. In this "moving grid" method, the fluid grid motion is driven by the motion of its boundaries, which are typically common boundaries of the moving fluid domain and the moving structure. This method introduces a fluid grid velocity *<sup>g</sup> v* <sup>G</sup> into the flow equations. When integrating Equation (2) over a fluid volume *V*, of which the surface *S* is moving with grid velocity *<sup>g</sup> v* <sup>G</sup> and has outward normal *<sup>n</sup>* <sup>G</sup> , Equation (2) becomes:

$$\frac{\partial}{\partial t}\int\_{V}dV + \int\_{S}\left(\vec{v} - \vec{v}\_{\mathcal{g}}\right) \cdot \vec{n}\,dS = 0\tag{7a}$$

$$\frac{\partial}{\partial t} \int\_{V} \rho \vec{v} \, dV + \int\_{S} \rho \vec{v} \left( \vec{v} - \vec{v}\_{\mathcal{S}} \right) \cdot \vec{n} \, dS = - \int\_{S} p \vec{n} \, dS + \int\_{S} \mu \left( \nabla \vec{v} + \left( \nabla \vec{v} \right)^{T} \right) \cdot \vec{n} \, dS \tag{7b}$$

When *<sup>g</sup> v v* <sup>=</sup> G G , this results in a purely Lagrangian description. When 0 *<sup>g</sup> <sup>v</sup>* <sup>=</sup> <sup>G</sup> <sup>G</sup> , a purely Eulerian description is recovered. Moreover, the fluid grid velocity is called "arbitrary" because it does not have to correspond to the fluid velocity. However, when the grid deformation becomes too large, as is the case with BMHVs, this could deteriorate the grid quality. Therefore, local remeshing is needed between time-steps. The main advantage of the ALE approach is its accuracy, because the grid is aligned with the fluid-structure interface. However, the use of remeshing (and thus interpolation) introduces artificial diffusivity and can become expensive for complex three-dimensional geometries. Several studies have used the ALE approach to simulate the dynamics of the ATS Open PivotTM Standard Heart Valve (Dumont et al., 2005, 2007), the St. Jude MedicalTM BMHV (Penrose et al., 2002; Redaelli et al., 2004; Dumont et al., 2007; Guivier et al., 2007, 2009; Nobili et al., 2007, 2008; Choi et al., 2009; Hong et al., 2009), and other valve types (Makhijani et al., 1997; Vierendeels et al., 2005, 2007; Bang et al., 2006; Morsi et al., 2007).

#### **2.2 Monolithic solver versus partitioned solver**

32 Aortic Valve

In the remainder of this section, classifications of FSI methods used in literature to simulate a BMHV will be made, and the characteristics of each of the mentioned methods will be

A first classification concerns the kinematical description of the domain. For the structural problem, the motion is usually described by the Lagrangian method, where the computational grid moves with the material velocity. This is in contrast to a fluid domain, in which the motion is generally described by the Eulerian method and, therefore, the computational grid does not deform. In case of FSI, both methods can be combined in

One approach is the "fixed grid" method, in which an immersed structure is allowed to "fictitiously" move through the Eulerian fluid grid in a Lagrangian way. The influence of the structure boundary on the fluid outside the structure is calculated by introducting body force sources into the Navier-Stokes equations, while keeping the velocity of the fictitiously overlapped fluid coupled to the structural velocity. Since the fluid grid is fixed, there is no need for remeshing and grid adaption. However, because the fluid-structure interface is not necessarily aligned with the spatial discretization and the data are, as such, interpolated, the flow field (and thus shear stresses) at this interface is not accurately calculated. Several fixed grid methods have been developed, such as the immersed boundary (IB) method, first proposed by Peskin (1972) for the simulation of heart valves. Borazjani et al. (2008) used the sharp interface CURVIB-method for simulating a BMHV. Other IB techniques were developed and used by Tai et al. (2007), De Tullio et al. (2009) and Xia et al. (2009). Also, the fictitious domain (FD) method can be used to simulate flexible heart valves (De Hart et al., 2000, 2003; Diniz dos Santos et al., 2008; Astorino et al., 2009). This fixed grid method uses Lagrange multipliers to impose the kinematic constraints. Van Loon et al. (2004) improved the accuracy of the FD method at the fluid-structure interface by proposing a local fluid grid

Another approach is to use the arbitrary Lagrangian-Eulerian (ALE) method for the fluid domain. In this "moving grid" method, the fluid grid motion is driven by the motion of its boundaries, which are typically common boundaries of the moving fluid domain and the

When integrating Equation (2) over a fluid volume *V*, of which the surface *S* is moving with

( ) <sup>0</sup> *<sup>g</sup> V S dV v v n dS*

*<sup>g</sup> V S S S v dV v v v n dS pn dS v v n dS <sup>t</sup>*

Eulerian description is recovered. Moreover, the fluid grid velocity is called "arbitrary" because it does not have to correspond to the fluid velocity. However, when the grid deformation becomes too large, as is the case with BMHVs, this could deteriorate the grid quality. Therefore, local remeshing is needed between time-steps. The main advantage of

<sup>∂</sup> + − ⋅ =− + ∇ + ∇ ⋅

When *<sup>g</sup> v v* <sup>=</sup> G G , this results in a purely Lagrangian description. When 0 *<sup>g</sup> <sup>v</sup>* <sup>=</sup>

( ) ( ( ) ) *<sup>T</sup>*

<sup>∂</sup> ∫ ∫ ∫∫ <sup>G</sup> GG G G G G G G (7b)

<sup>G</sup> , Equation (2) becomes:

<sup>+</sup> −⋅ = <sup>∂</sup> ∫ ∫ <sup>G</sup> G G (7a)

μ

<sup>G</sup> into the flow equations.

<sup>G</sup> <sup>G</sup> , a purely

**2.1 Fixed grid techniques versus moving grid techniques** 

several approaches in order to describe the motion of the domain.

moving structure. This method introduces a fluid grid velocity *<sup>g</sup> v*

*t* ∂

<sup>G</sup> and has outward normal *<sup>n</sup>*

 ρ

adaption at the structure boundary.

ρ

grid velocity *<sup>g</sup> v*

explained.

Secondly, one can classify each FSI simulation by using a partitioned solver or a monolithic solver. In the monolithic approach, the entire FSI problem is simultaneously simulated by one solver.

This is in contrast to the partitioned approach, which solves the flow and the structural problem separately and, therefore, mostly uses different specialized solvers. The partitioned approach is used to simulate heart valves in Makhijani et al. (1997), Penrose et al. (2002), Redaelli et al. (2004), Dumont et al. (2005, 2007), Vierendeels et al. (2005, 2007), Bang et al. (2006), Guivier et al. (2007, 2009), Nobili et al. (2007, 2008), Morsi et al. (2007), Tai et al. (2007), Borazjani et al. (2008), Diniz dos Santos et al. (2008), Astorino et al. (2009), Choi et al. (2009), De Tullio et al. (2009), Hong et al. (2009) and, finally, Xia et al. (2009). In order to obtain the interaction between the fluid and the structure, data exchange at the fluid-structure interface and a coupling scheme between the separated solvers are needed. Unfortunately, not every coupling scheme converges quickly. The instability of coupling schemes without relaxation is analytically explained in Vierendeels et al. (2005) and Borazjani et al. (2008) for the case of BMHVs, and it is also demonstrated by the flow through arteries (Causin et al., 2005; Degroote et al., 2008, 2010). It is concluded that relaxation can be used to obtain stable and efficient approximations for the subsequent angular accelerations of the valve leaflets. Several coupling schemes with relaxation have thus been developed, and they can be divided into loose and strong coupling schemes.

In the loose coupling methods, only one coupling iteration is needed in each time-step, since the solution of the flow field at time-step *n* is used to calculate the angular accelerations of the leaflets for the next time-step *n*+1:

$$
\dot{\theta}\_{i}^{n+1} = \left(1 - \alpha\_{i}^{n+1}\right) \cdot \dot{\theta}\_{i}^{n} + \alpha\_{i}^{n+1} \cdot \frac{M\_{i}^{n}}{I\_{i}} \tag{8}
$$

A relaxation factor ( *<sup>n</sup>* <sup>1</sup> *i* ω <sup>+</sup> < 1) is necessary, since the scheme without relaxation ( *<sup>n</sup>* <sup>1</sup> *i* ω <sup>+</sup> = 1) is unstable (Causin et al., 2005; Vierendeels et al., 2005; Borazjani et al., 2008), as mentioned above. The loose coupling formulation is often used to simulate heart valves (Redaelli et al., 2004; Morsi et al., 2007; Nobili et al., 2007; Tai et al., 2007; Borazjani et al., 2008; Xia et al., 2009). It has the main benefit of a low computational cost because of the lack of a coupling iteration loop within each time-step. However, this lack implies that dynamic equilibrium at the fluid-structure interface (Equation (6)) is not necessarily achieved, which leads to unphysical oscillations in the leaflet movement and the flow and pressure field, as described in Annerel et al. (2011).

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 35

( ) <sup>ˆ</sup> 1, 1 1, 1, 1, 1, <sup>1</sup> *nk nk nk nk nk M IM i ii*

1, 1 1, 1, 1, 1

*nk nk nk nk i i ii i*

To simulate a BMHV in a partitioned way, several methods can be used to calculate an

For the fixed relaxation factor, the factor value is kept constant during the entire simulation,

 ω

Fixed relaxation was used to simulate the dynamics of a BMHV with a partitioned solver by Penrose et al. (2002), Nobili et al. (2007, 2008), Borazjani et al. (2008), De Tullio et al. (2009), and Hong et al. (2009). The main disadvantage of such a fixed relaxation is the lack of a physical meaning of the relaxation factor (Annerel et al., 2011). Therefore, the selection of an appropriate factor value is ad hoc and will be done primarily through trial-and-error, as

For the dynamic relaxation factor, the factor value is updated in each coupling iteration of each time-step. Typically, the Aitken Δ2 relaxation is used, as described in Irons et al. (1969), Mok et al. (2001), and Küttler et al. (2008). The Aitken ∆2 relaxation updates the value of the

*<sup>i</sup> <sup>T</sup> nk nk nk nk*

− − = = −

( ) ( )

+ +− + +−

*<sup>T</sup> nk nk nk nk*

*θ θ r r*

*rr rr*

1, 1, 1 1, 1, 1

1, 1, 1 1, 1, 1

1, 1 1, 1

*n k*

θ

<sup>+</sup> <sup>+</sup>

<sup>+</sup> <sup>+</sup> <sup>⎡</sup> <sup>⎤</sup> <sup>⎢</sup> <sup>−</sup> <sup>⎥</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>=</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>⎢</sup> <sup>−</sup> <sup>⎥</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>⎣</sup> <sup>⎦</sup>

θ

*n k*

*r* (15)

*n k*

1, 2 1, 2

*n k*

2

( ) ( )

*M I M I*

− −

1, 1

*n k*

+

Partitioned simulations of heart valves using the Aitken Δ2 relaxation method are reported

More recently, however, a (quasi-Newton) method with a dynamically changing relaxation matrix has been developed (Annerel et al., 2010; Dahl et al., 2010) and subsequently optimized in Annerel et al. (2011). The method is an extension of Vierendeels et al. (2005)

through a linearization of the dynamic equilibrium. Thus, while taking into account the mutual interaction between the leaflets, Equation (6) is linearized for each coupling iteration

in Borazjani et al. (2008), Diniz dos Santos et al. (2008), and Astorino et al. (2009).

and predicts the moments of the next coupling iteration (i.e. <sup>ˆ</sup> *n k* 1, 1 *Mi*

+ +− + +−

θω *i ii* ++ + + + + =− ⋅ + ⋅ (11)

*n k*

*i*

*<sup>i</sup> Cst* <sup>+</sup> = = (13)

(14)

<sup>+</sup> <sup>+</sup> in Equation (10))

*M I*

<sup>+</sup> that stabilizes the solution process. In the following,

<sup>+</sup> ++ + + + =− ⋅ + ⋅ (12)

( ) 1,

 ω

ω

 ωθ

using a fixed relaxation factor and a dynamic relaxation factor is discussed.

as described in Le Tallec & Mouro (2001) for an industrial shock absorber valve:

*n k* 1, ω

*i* ω

Inserted in Equation (10):

θ

factor in each coupling iteration of each time-step *n*+1:

1, 1,

 ω

1, 1

*n k*

+

1

*n ,k*

+

1 2

*n ,k θ θ*

+ <sup>⎡</sup> <sup>⎤</sup> <sup>=</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>⎢</sup> <sup>⎥</sup> <sup>⎣</sup> <sup>⎦</sup> *<sup>θ</sup>* and

*nk nk*

+ +

ω

appropriate relaxation factor *n k* 1,

noted by De Tullio et al. (2009).

with

*k+*1 of time-step *n+*1:

Dynamic equilibrium at the fluid-structure interface can be obtained by introducing a coupling iteration loop within each time-step, as is the case in the strong coupling methods. Generally, each of the strong coupling iterations follows the same pattern, as visualized in Fig. 2. At the beginning of each coupling iteration *k* of time-step *n*+1, the motion of the leaflets is computed from the angular accelerations *n k* 1, θ*i* <sup>+</sup> . The mesh is moved and the flow equations are solved. From the flow field, the moments *Mi <sup>n</sup>*+1,*<sup>k</sup>* are calculated. Finally, the convergence of the dynamic equilibrium at the fluid-structure interface, expressed by Equation (6), is checked. When this dynamic equilibrium is obtained, a new time-step is initiated. However, when dynamic equilibrium is not achieved, a new coupling iteration *k*+1 is initiated, and thus new angular accelerations *n k* 1, 1 *i* θ <sup>+</sup> <sup>+</sup> need to be calculated. Therefore, introducing a coupling iteration loop requires, in each coupling iteration *k* of time-step *n*+1, a stable and efficient approximation of the angular accelerations for the next coupling iteration *k*+1.

Fig. 2. Simplified flow chart of a strong FSI coupling algorithm with two degrees of freedom. *n* = time-step, *k* = coupling iteration step, *i* = leaflet number

A strong coupling scheme without relaxation can easily be proposed. From the moments of coupling iteration *k* in time-step *n*+1, the angular accelerations of the next coupling iteration *k*+1, i.e. *n k* 1, 1 θ*i* + + , can be calculated:

$$
\dot{\theta}\_i^{n+1,k+1} = \frac{M\_i^{n+1,k}}{I\_i} \tag{9}
$$

Unfortunately, such fixed-point iterations, also called Gauss-Seidel iterations without relaxation, are unstable for BMHVs (Vierendeels et al., 2005). Therefore, the scheme needs to be stabilized by using an appropriate prediction of the moments of the next coupling iteration *k*+1:

$$
\ddot{\theta}\_{i}^{n+1,k+1} = \frac{\hat{M}\_{i}^{n+1,k+1}}{I\_{i}} \tag{10}
$$

with <sup>ˆ</sup> *n k* 1, 1 *Mi* + + denoting the predicted moment. Several methods are used in the literature to calculate a stable value for <sup>ˆ</sup> *n k* 1, 1 *Mi* <sup>+</sup> <sup>+</sup> . Usually, this is achieved using a relaxation scheme, which leads to fixed-point iterations with relaxation ( *n k* 1, *i* ω <sup>+</sup> < 1), also called the Gauss-Seidel coupling method with relaxation:

$$
\hat{M}\_i^{n+1,k+1} = \left(1 - \alpha\_i^{n+1,k}\right) \cdot I\_i \ddot{\theta}\_i^{n+1,k} + \alpha\_i^{n+1,k} \cdot M\_i^{n+1,k} \tag{11}
$$

Inserted in Equation (10):

34 Aortic Valve

Dynamic equilibrium at the fluid-structure interface can be obtained by introducing a coupling iteration loop within each time-step, as is the case in the strong coupling methods. Generally, each of the strong coupling iterations follows the same pattern, as visualized in Fig. 2. At the beginning of each coupling iteration *k* of time-step *n*+1, the motion of the

convergence of the dynamic equilibrium at the fluid-structure interface, expressed by Equation (6), is checked. When this dynamic equilibrium is obtained, a new time-step is initiated. However, when dynamic equilibrium is not achieved, a new coupling iteration *k*+1

introducing a coupling iteration loop requires, in each coupling iteration *k* of time-step *n*+1, a stable and efficient approximation of the angular accelerations for the next coupling

Solve flow equations

Dyn. equilibrium?

Fig. 2. Simplified flow chart of a strong FSI coupling algorithm with two degrees of freedom.

A strong coupling scheme without relaxation can easily be proposed. From the moments of coupling iteration *k* in time-step *n*+1, the angular accelerations of the next coupling iteration

1, 1

*i*

θ

*i*

θ

which leads to fixed-point iterations with relaxation ( *n k* 1,

*n k i*

Unfortunately, such fixed-point iterations, also called Gauss-Seidel iterations without relaxation, are unstable for BMHVs (Vierendeels et al., 2005). Therefore, the scheme needs to be stabilized by using an appropriate prediction of the moments of the next coupling

> 1, 1 <sup>ˆ</sup> *n k n k i*

*M I*

*kni <sup>n</sup>*

*<sup>n</sup>*+1,*<sup>k</sup>* from flow field

1,

1, 1

*i*

+ + denoting the predicted moment. Several methods are used in the literature to

<sup>+</sup> + + <sup>=</sup> (9)

<sup>+</sup> <sup>+</sup> + + <sup>=</sup> (10)

<sup>+</sup> < 1), also called the Gauss-

<sup>+</sup> <sup>+</sup> . Usually, this is achieved using a relaxation scheme,

*i* ω

*n k*

*i*

*M I*

Calculate new time *tn+1* = *tn*+Δ*tn+1* and estimate *k*= 0

*kn i* + ,1 θ

Calculate *Mi*

[no] [yes] *kn* ++ 1,1 *<sup>i</sup>*

θ*i*

*i* θ <sup>+</sup> . The mesh is moved and the flow

<sup>+</sup> <sup>+</sup> need to be calculated. Therefore,

*n*+ 0,1 θ*i* 

*<sup>n</sup>*+1,*<sup>k</sup>* are calculated. Finally, the

*n*= *n*+1

*<sup>i</sup>* <sup>1</sup> <sup>+</sup> ,1 <sup>=</sup> <sup>+</sup> θ

θ

leaflets is computed from the angular accelerations *n k* 1,

is initiated, and thus new angular accelerations *n k* 1, 1

*k* = *k*+1 Mesh movement with

*n* = time-step, *k* = coupling iteration step, *i* = leaflet number

iteration *k*+1.

Calculate

*k*+1, i.e. *n k* 1, 1 θ*i*

iteration *k*+1:

with <sup>ˆ</sup> *n k* 1, 1 *Mi*

θ

+ + , can be calculated:

calculate a stable value for <sup>ˆ</sup> *n k* 1, 1 *Mi*

Seidel coupling method with relaxation:

equations are solved. From the flow field, the moments *Mi*

$$
\ddot{\theta}\_{i}^{n+1,k+1} = \left(1 - \alpha\_{i}^{n+1,k}\right) \cdot \ddot{\theta}\_{i}^{n+1,k} + \alpha\_{i}^{n+1,k} \cdot \frac{M\_{i}^{n+1,k}}{I\_{i}} \tag{12}
$$

To simulate a BMHV in a partitioned way, several methods can be used to calculate an appropriate relaxation factor *n k* 1, *i* ω <sup>+</sup> that stabilizes the solution process. In the following, using a fixed relaxation factor and a dynamic relaxation factor is discussed.

For the fixed relaxation factor, the factor value is kept constant during the entire simulation, as described in Le Tallec & Mouro (2001) for an industrial shock absorber valve:

$$
\alpha\_i^{n+1,k} = \alpha = \text{Cst} \tag{13}
$$

Fixed relaxation was used to simulate the dynamics of a BMHV with a partitioned solver by Penrose et al. (2002), Nobili et al. (2007, 2008), Borazjani et al. (2008), De Tullio et al. (2009), and Hong et al. (2009). The main disadvantage of such a fixed relaxation is the lack of a physical meaning of the relaxation factor (Annerel et al., 2011). Therefore, the selection of an appropriate factor value is ad hoc and will be done primarily through trial-and-error, as noted by De Tullio et al. (2009).

For the dynamic relaxation factor, the factor value is updated in each coupling iteration of each time-step. Typically, the Aitken Δ2 relaxation is used, as described in Irons et al. (1969), Mok et al. (2001), and Küttler et al. (2008). The Aitken ∆2 relaxation updates the value of the factor in each coupling iteration of each time-step *n*+1:

$$\alpha\_{i}^{n+1,k} = \alpha^{n+1,k} = -\frac{\left(\ddot{\theta}^{n+1,k} - \ddot{\theta}^{n+1,k-1}\right)^{T} \left(r^{n+1,k} - r^{n+1,k-1}\right)}{\left(r^{n+1,k} - r^{n+1,k-1}\right)^{T} \left(r^{n+1,k} - r^{n+1,k-1}\right)}\tag{14}$$

with

$$\ddot{\theta}^{n+1,k} = \begin{bmatrix} \ddot{\theta}\_1^{n+1,k} \\ \ddot{\theta}\_2^{n+1,k} \end{bmatrix} \text{ and } \quad \mathbf{r}^{n+1,k} = \begin{bmatrix} M\_1^{n+1,k} \\ I\_1 \\ \frac{M\_1^{n+1,k}}{I\_2} - \ddot{\theta}\_1^{n+1,k} \\ \frac{M\_2^{n+1,k}}{I\_2} - \ddot{\theta}\_2^{n+1,k} \end{bmatrix} \tag{15}$$

Partitioned simulations of heart valves using the Aitken Δ2 relaxation method are reported in Borazjani et al. (2008), Diniz dos Santos et al. (2008), and Astorino et al. (2009).

More recently, however, a (quasi-Newton) method with a dynamically changing relaxation matrix has been developed (Annerel et al., 2010; Dahl et al., 2010) and subsequently optimized in Annerel et al. (2011). The method is an extension of Vierendeels et al. (2005) and predicts the moments of the next coupling iteration (i.e. <sup>ˆ</sup> *n k* 1, 1 *Mi* <sup>+</sup> <sup>+</sup> in Equation (10)) through a linearization of the dynamic equilibrium. Thus, while taking into account the mutual interaction between the leaflets, Equation (6) is linearized for each coupling iteration *k+*1 of time-step *n+*1:

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 37

transition to a small-scale turbulent state downstream of the valve (Borazjani et al., 2008;

The valve leaflets start to close at the beginning of the steepest flow deceleration (De Tullio et al., 2009). Since the leaflets need to rotate over a large angle, some regurgitation occurs. The total volume of this reverse flow is denoted as the "closing volume" (as also described

During leaflet closure, the angular velocity of the leaflets keeps increasing until the closed position is reached. This gives rise to very large angular velocities (and thus stresses) when the leaflets impact the blocking mechanism at the fully closed position. Therefore, the closing kinematics are very different from those at the leaflet opening phase because the end of the opening occurs with small angular velocities, while at the end of the closing, the

During the closing phase, the decelerating flow field is governed by small-scale eddies and

The leaflets reach the fully closed position at the negative peak of the flow rate (De Tullio et al., 2009). Due to the negative transvalvular pressure gradient when the leaflets are closed, flow leaks through the small gaps between the leaflets and through the gaps between the leaflets and valve casing (in particular, near the hinges), giving rise to squeezed jet flow. The total amount of this regurgitant flow is denoted as the "leakage volume" (see Section 5). After valve closure, the turbulent structures in the flow slowly decay. Subsequently, the residual vortical structures are rapidly washed out at the beginning of a new cycle by the incoming accelerated flow when the valve reopens. (Borazjani et al., 2008; De Tullio et al.,

**4. Three-dimensional strongly coupled partitioned FSI simulations of a BMHV**  To illustrate the previous section, a BMHV is simulated in three different geometries. The used BMHV is a simplified model of the 25-mm ATS Open PivotTM Standard Heart Valve in aortic position, with the orifice inner diameter measuring 20.8mm. The valve is simplified at the hinge regions by cutting away the blocking mechanism and hinges at the casing. Because of this simplification, the resulting opening velocity of the valve leaflets could become slightly overestimated since the additional counteracting moment created by the decelerated squeeze flows near the pivot hinge regions (i.e. the so-called pivot effect, as

The valve is subsequently placed in three geometries. The first geometry consists of a rigid straight tube, as visualized in Fig. 3(*left*). The second geometry also consists of a rigid straight tube upstream of the valve, but rigid sinuses of Valsalva are added downstream of the valve. Such sinuses of Valsalva are anatomically present in the ascending aortic root and influence the valve closing (Grigioni et al., 2004). The sinuses of Valsalva are based on the geometry described in Reul et al. (1990) and are positioned symmetrically with respect to the leaflet rotation axes (Fig. 3(*middle*)). In the third geometry, the same sinuses of Valsalva are used, but they are placed asymmetrically (angle: 30°) with respect to the leaflet rotation

axes in such a way that one of the leaflets directly faces one sinus (Fig. 3(*right*)).

angular velocity attains peak value (De Tullio et al., 2009).

observed in the experiments in Feng et al. (1999)) is absent.

turbulent vortices (Borazjani et al., 2008; Sotiropoulos et al., 2009).

Sotiropoulos et al., 2009).

**Closing phase** 

in Section 5).

**Fully closed position** 

2009; Sotiropoulos et al., 2009).

$$\begin{cases} M\_1^{n+1,k} + \left(\frac{\partial M\_1}{\partial \ddot{\theta}\_1}\right)^{n+1,k} \left(\ddot{\theta}\_1^{n+1,k+1} - \ddot{\theta}\_1^{n+1,k}\right) + \left(\frac{\partial M\_1}{\partial \ddot{\theta}\_2}\right)^{n+1,k} \left(\ddot{\theta}\_2^{n+1,k+1} - \ddot{\theta}\_2^{n+1,k}\right) = I\_1 \cdot \ddot{\theta}\_1^{n+1,k+1} \\\ M\_2^{n+1,k} + \left(\frac{\partial M\_2}{\partial \ddot{\theta}\_1}\right)^{n+1,k} \left(\ddot{\theta}\_1^{n+1,k+1} - \ddot{\theta}\_1^{n+1,k}\right) + \left(\frac{\partial M\_2}{\partial \ddot{\theta}\_2}\right)^{n+1,k} \left(\ddot{\theta}\_2^{n+1,k+1} - \ddot{\theta}\_2^{n+1,k}\right) = I\_2 \cdot \ddot{\theta}\_2^{n+1,k+1} \end{cases} \tag{16}$$

The components of the Jacobian are the derivatives of the moments (exerted by the flow on the leaflets) with respect to changes in leaflet angular accelerations. This Jacobian is approximated with finite differences and is numerically calculated from the flow solver by variations of the leaflet angular accelerations. This method outperforms the fixed relaxation and Aitken Δ2 relaxation in needed coupling iterations per time-step and CPU time (Annerel et al., 2010, 2011).

## **3. Insights into the general dynamics and flow fields of an aortic BMHV**

The dynamics of a BMHV depend on passive movement. Therefore, the opening and closure of the leaflets are governed by the pressure gradients and flow fields in the heart and arteries (in case of atrioventricular valves) (Butany et al., 2003).

In the following, the kinematics and dynamics of a BMHV are discussed and obtained by numerical simulations in an axisymmetric geometry (Yoganathan et al., 2004; Borazjani et al., 2008; De Tullio et al., 2009; Sotiropoulos et al., 2009) and verified by experiments (Dasi et al., 2007). These discussed dynamics are also illustrated in the numerical simulations done in Section 4.

#### **Opening phase**

The contraction of the left ventricle at the beginning of systole induces an increase of the left ventricular pressure. Because of the resulting positive transvalvular pressure gradient, the flow starts to accelerate and induces the opening of the valve.

In a BMHV, the valve leaflets initially open with a large increase in angular acceleration (and related angular velocity). However, as the valve opens, the leaflets tend to align with the axial flow, and their local linear velocity tends to become orthogonal to the main flow stream lines, which produces a pressure moment that decreases the angular acceleration and lowers the angular velocity (De Tullio et al., 2009). This deceleration is beneficial because it results in a very small angular velocity of the leaflets when the leaflets reach the fully open position; therefore, it significantly lowers the impact forces at the blocking mechanism.

During the early leaflet opening phase, three jet flows are formed, with the roll-up of the valve housing shear layer into the aortic sinuses of Valsalva and the formation of two shear layers shed from the tips of the valve leaflets (Borazjani et al., 2008; De Tullio et al., 2009; Sotiropoulos et al., 2009). Subsequently, these leaflet shear layers break down and largescale von Karman-like vortex shedding emerges (Sotiropoulos et al., 2009).

#### **Fully open position**

When the leaflets have reached the fully open position, the flow rate continues increasing until its maximum value at peak systole. Due to the acceleration of the flow, the formation of small-scale turbulence is prevented, and the bulk of the flow remains laminar (De Tullio et al., 2009). However, at peak systole, the flow starts to decelerate. Subsequently, the largescale vortices in the sinuses of Valsalva and the leaflet shear layers rapidly undergo the transition to a small-scale turbulent state downstream of the valve (Borazjani et al., 2008; Sotiropoulos et al., 2009).

#### **Closing phase**

36 Aortic Valve

*n k nk nk nk nk nk*

+ − ⎜ ⎟ + − ⎜ ⎟ <sup>=</sup> <sup>⋅</sup> ∂ ∂ ⎝ ⎠ ⎝ ⎠

+ − ⎜ ⎟ + − ⎜ ⎟ <sup>=</sup> <sup>⋅</sup> ∂ ∂ ⎝ ⎠ ⎝ ⎠

The components of the Jacobian are the derivatives of the moments (exerted by the flow on the leaflets) with respect to changes in leaflet angular accelerations. This Jacobian is approximated with finite differences and is numerically calculated from the flow solver by variations of the leaflet angular accelerations. This method outperforms the fixed relaxation and Aitken Δ2 relaxation in needed coupling iterations per time-step and CPU time (Annerel

1, 1,

*M M M I*

*M M M I*

**3. Insights into the general dynamics and flow fields of an aortic BMHV** 

The dynamics of a BMHV depend on passive movement. Therefore, the opening and closure of the leaflets are governed by the pressure gradients and flow fields in the heart and

In the following, the kinematics and dynamics of a BMHV are discussed and obtained by numerical simulations in an axisymmetric geometry (Yoganathan et al., 2004; Borazjani et al., 2008; De Tullio et al., 2009; Sotiropoulos et al., 2009) and verified by experiments (Dasi et al., 2007). These discussed dynamics are also illustrated in the numerical simulations done

The contraction of the left ventricle at the beginning of systole induces an increase of the left ventricular pressure. Because of the resulting positive transvalvular pressure gradient, the

In a BMHV, the valve leaflets initially open with a large increase in angular acceleration (and related angular velocity). However, as the valve opens, the leaflets tend to align with the axial flow, and their local linear velocity tends to become orthogonal to the main flow stream lines, which produces a pressure moment that decreases the angular acceleration and lowers the angular velocity (De Tullio et al., 2009). This deceleration is beneficial because it results in a very small angular velocity of the leaflets when the leaflets reach the fully open position; therefore, it significantly lowers the impact forces at the blocking mechanism. During the early leaflet opening phase, three jet flows are formed, with the roll-up of the valve housing shear layer into the aortic sinuses of Valsalva and the formation of two shear layers shed from the tips of the valve leaflets (Borazjani et al., 2008; De Tullio et al., 2009; Sotiropoulos et al., 2009). Subsequently, these leaflet shear layers break down and large-

When the leaflets have reached the fully open position, the flow rate continues increasing until its maximum value at peak systole. Due to the acceleration of the flow, the formation of small-scale turbulence is prevented, and the bulk of the flow remains laminar (De Tullio et al., 2009). However, at peak systole, the flow starts to decelerate. Subsequently, the largescale vortices in the sinuses of Valsalva and the leaflet shear layers rapidly undergo the

*n k n k*

+ +

1, 1, 1, 2 2 1, 1 1, 1, 1 1, 2 1 1 2 2 2

*n k n k n k nk nk nk nk*

+ + + ++ + ++ +

1 2

θθ

⎛ ⎞ ∂ ∂⎛ ⎞

⎛ ⎞ ∂ ∂⎛ ⎞

arteries (in case of atrioventricular valves) (Butany et al., 2003).

flow starts to accelerate and induces the opening of the valve.

scale von Karman-like vortex shedding emerges (Sotiropoulos et al., 2009).

θ

θ

⎧ ⎪ ⎪⎪ ⎨ ⎪ ⎪ ⎪⎩

et al., 2010, 2011).

in Section 4. **Opening phase** 

**Fully open position** 

θ

1 2

 θ

( ) ( )

θθ

θθ

 θ

> θ

2 *n k* + + (16)

1, 1 1 1, 1 1, 1, 1 1, 1, 1 1 1 1 2 2 1 1

+ ++ + ++ + + +

 θ

> θ

( ) ( )

1, 1

The valve leaflets start to close at the beginning of the steepest flow deceleration (De Tullio et al., 2009). Since the leaflets need to rotate over a large angle, some regurgitation occurs. The total volume of this reverse flow is denoted as the "closing volume" (as also described in Section 5).

During leaflet closure, the angular velocity of the leaflets keeps increasing until the closed position is reached. This gives rise to very large angular velocities (and thus stresses) when the leaflets impact the blocking mechanism at the fully closed position. Therefore, the closing kinematics are very different from those at the leaflet opening phase because the end of the opening occurs with small angular velocities, while at the end of the closing, the angular velocity attains peak value (De Tullio et al., 2009).

During the closing phase, the decelerating flow field is governed by small-scale eddies and turbulent vortices (Borazjani et al., 2008; Sotiropoulos et al., 2009).

#### **Fully closed position**

The leaflets reach the fully closed position at the negative peak of the flow rate (De Tullio et al., 2009). Due to the negative transvalvular pressure gradient when the leaflets are closed, flow leaks through the small gaps between the leaflets and through the gaps between the leaflets and valve casing (in particular, near the hinges), giving rise to squeezed jet flow. The total amount of this regurgitant flow is denoted as the "leakage volume" (see Section 5).

After valve closure, the turbulent structures in the flow slowly decay. Subsequently, the residual vortical structures are rapidly washed out at the beginning of a new cycle by the incoming accelerated flow when the valve reopens. (Borazjani et al., 2008; De Tullio et al., 2009; Sotiropoulos et al., 2009).

#### **4. Three-dimensional strongly coupled partitioned FSI simulations of a BMHV**

To illustrate the previous section, a BMHV is simulated in three different geometries. The used BMHV is a simplified model of the 25-mm ATS Open PivotTM Standard Heart Valve in aortic position, with the orifice inner diameter measuring 20.8mm. The valve is simplified at the hinge regions by cutting away the blocking mechanism and hinges at the casing. Because of this simplification, the resulting opening velocity of the valve leaflets could become slightly overestimated since the additional counteracting moment created by the decelerated squeeze flows near the pivot hinge regions (i.e. the so-called pivot effect, as observed in the experiments in Feng et al. (1999)) is absent.

The valve is subsequently placed in three geometries. The first geometry consists of a rigid straight tube, as visualized in Fig. 3(*left*). The second geometry also consists of a rigid straight tube upstream of the valve, but rigid sinuses of Valsalva are added downstream of the valve. Such sinuses of Valsalva are anatomically present in the ascending aortic root and influence the valve closing (Grigioni et al., 2004). The sinuses of Valsalva are based on the geometry described in Reul et al. (1990) and are positioned symmetrically with respect to the leaflet rotation axes (Fig. 3(*middle*)). In the third geometry, the same sinuses of Valsalva are used, but they are placed asymmetrically (angle: 30°) with respect to the leaflet rotation axes in such a way that one of the leaflets directly faces one sinus (Fig. 3(*right*)).

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 39

(a)

(b)

(c) Fig. 4. Angular position of the leaflets (relative to maximal opening angle) and the aortic flow pulse velocity (*a*), angular velocity of the leaflets (*b*) and angular acceleration of the leaflets (*c*). The two leaflets of the geometry with the straight tube perform a complete symmetrical movement. The time levels at which the contours in Fig. 5 are shown are

visualized by the vertical lines

For the geometries, the upstream tube has a diameter of 22mm and is 75mm long. The downstream geometry is 95mm long. The diameter of the downstream tube is 27.36 mm for the sinuses of Valsalva and 22mm for the straight tube.

These geometries are based on clinical practice, because when implanting the BMHV, the surgeon can choose to preserve the sinuses of Valsalva or to cut them away and replace the entire ascending aortic root (in the so-called Bentall procedure (Bentall et al., 1968)). Moreover, the surgeon can choose to position the valve symmetrically to the physiological sinuses of Valsalva, or to position it asymmetrically.

An unstructured fluid grid, consisting of approximately 800 000 tetrahedral cells, is generated in the geometries. Two cell layers are generated in the gap (which measures 0.1mm) between the leaflets and the casing near the hinge region. The ALE approach is followed, which means that the fluid grid follows the motion of the structure and subsequently needs an update to maintain good mesh quality. This update is done using a remeshing method and spring-based smoothing.

Fig. 3. Isometric view of the different geometries. *Left*: first geometry with straight tube. *Middle*: second geometry with symmetrically placed sinuses of Valsalva. *Right*: third geometry with asymmetrically positioned sinuses of Valsalva downstream of the valve

The dynamics of the BMHV are modeled by the strongly coupled partitioned quasi-Newton algorithm that was recently developed (Annerel et al., 2011). An inlet aortic flow pulse with a period of 1s (displayed in Fig. 4(*a*)) is imposed upstream and was previously used in Dumont et al. (2005, 2007) and Annerel et al. (2010, 2011). The flow pulse profile is uniform. A physiological pressure profile is imposed at the downstream outlet boundary. Note, however, that in a rigid geometry the pressure level does not affect the flow field since only the pressure gradient appears in the equations.

Blood is modeled as an incompressible Newtonian fluid with density and viscosity respectively equal to 1050kg/m3 and 4E-3Pa·s (i.e. the high shear rate limit viscosity of blood). Although real blood is a heterogeneous non-Newtonian fluid, the modeling of blood as a homogeneous Newtonian fluid for high shear rates is widely agreed upon for flow in large arteries and valves (Paul et al., 2003; Sotiropoulos et al., 2009). Nevertheless, it is important to keep in mind that when studying low levels of shear stresses, for example in the recirculation regions and vortices in the wake, the non-Newtonian effects could become important and should be taken into account when assessing the hemodynamics of the valve (Sotiropoulos et al., 2009). In such cases it is valuable to model the blood as a non-Newtonian fluid as in the Cross and Carreau model (Cross, 1965; Carreau, 1972). However, modeling of the fine-flow features and the hemodynamics is beyond the scope of the described simulations and, therefore, a Newtonian blood model is used.

For the geometries, the upstream tube has a diameter of 22mm and is 75mm long. The downstream geometry is 95mm long. The diameter of the downstream tube is 27.36 mm for

These geometries are based on clinical practice, because when implanting the BMHV, the surgeon can choose to preserve the sinuses of Valsalva or to cut them away and replace the entire ascending aortic root (in the so-called Bentall procedure (Bentall et al., 1968)). Moreover, the surgeon can choose to position the valve symmetrically to the physiological

An unstructured fluid grid, consisting of approximately 800 000 tetrahedral cells, is generated in the geometries. Two cell layers are generated in the gap (which measures 0.1mm) between the leaflets and the casing near the hinge region. The ALE approach is followed, which means that the fluid grid follows the motion of the structure and subsequently needs an update to maintain good mesh quality. This update is done using a

Fig. 3. Isometric view of the different geometries. *Left*: first geometry with straight tube. *Middle*: second geometry with symmetrically placed sinuses of Valsalva. *Right*: third geometry with asymmetrically positioned sinuses of Valsalva downstream of the valve

The dynamics of the BMHV are modeled by the strongly coupled partitioned quasi-Newton algorithm that was recently developed (Annerel et al., 2011). An inlet aortic flow pulse with a period of 1s (displayed in Fig. 4(*a*)) is imposed upstream and was previously used in Dumont et al. (2005, 2007) and Annerel et al. (2010, 2011). The flow pulse profile is uniform. A physiological pressure profile is imposed at the downstream outlet boundary. Note, however, that in a rigid geometry the pressure level does not affect the flow field since only

Blood is modeled as an incompressible Newtonian fluid with density and viscosity respectively equal to 1050kg/m3 and 4E-3Pa·s (i.e. the high shear rate limit viscosity of blood). Although real blood is a heterogeneous non-Newtonian fluid, the modeling of blood as a homogeneous Newtonian fluid for high shear rates is widely agreed upon for flow in large arteries and valves (Paul et al., 2003; Sotiropoulos et al., 2009). Nevertheless, it is important to keep in mind that when studying low levels of shear stresses, for example in the recirculation regions and vortices in the wake, the non-Newtonian effects could become important and should be taken into account when assessing the hemodynamics of the valve (Sotiropoulos et al., 2009). In such cases it is valuable to model the blood as a non-Newtonian fluid as in the Cross and Carreau model (Cross, 1965; Carreau, 1972). However, modeling of the fine-flow features and the hemodynamics is beyond the scope of the

described simulations and, therefore, a Newtonian blood model is used.

the sinuses of Valsalva and 22mm for the straight tube.

sinuses of Valsalva, or to position it asymmetrically.

remeshing method and spring-based smoothing.

the pressure gradient appears in the equations.

Fig. 4. Angular position of the leaflets (relative to maximal opening angle) and the aortic flow pulse velocity (*a*), angular velocity of the leaflets (*b*) and angular acceleration of the leaflets (*c*). The two leaflets of the geometry with the straight tube perform a complete symmetrical movement. The time levels at which the contours in Fig. 5 are shown are visualized by the vertical lines

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 41

Since the regurgitation lowers the net cardiac output and therefore enlarges the workload of

(a)

(b)

(c)

(d)

perpendicular to the leaflet rotation axes, for the geometry with the straight tube (*left*), the symmetrically placed sinuses of Valsalva (*middle*) and the asymmetrical sinuses of Valsalva (*right*). The plots are taken at *t* = 0.025s (*a*), *t* = 0.125s (*b*), *t* = 0.250s (*c*), and *t* = 0.375s (*d*),

Fig. 5. Velocity Magnitude Contours in m/s visualized on a longitudinal section

represented by the vertical lines in Fig. 4

the heart, its percentage should be minimized when designing heart valves.

No turbulence model is used, thus implying laminar flow. A no-slip boundary condition is applied at the walls. The valve is initially set in the closed position. The moment of inertia of one rigid valve leaflet about its rotation axis is equal to 9.94E-9kg·m2.

The results of the simulations are depicted in Fig. 4 and Fig. 5. The angular positions of the leaflets are presented in Fig. 4(*a*), relative to the fully opened position. Therefore, 0 and 1 refer, respectively, to the fully closed and fully open position.

Although the valve leaflets open completely in the first geometry with the straight tube, the results show that this maximum opening position is not reached in both geometries with the sinuses of Valsalva. Such incomplete opening of the ATS Open PivotTM Standard Heart Valve in a divergent geometry is explained by the greater sensitivity of the leaflet movement to the flow field compared with other BMHV designs, since the leaflets extend farther in the flow downstream of the orifice than is the case in other valve designs (Feng et al., 1999). Therefore, the valve does not open completely in the divergent transvalvular flow caused by the enlargement of the sinuses of Valsalva because the leaflets tend to align with the streamlines. In the straight tube, however, the valve leaflets open completely.

This incomplete opening also explains the difference in the closing phase between the geometries, since the valve in the sinuses of Valsalva geometries is closed sooner. This is because the leaflets in the straight tube reach the completely open position and therefore need to rotate over a greater angle in order to close. Hence, they have a greater closing volume (Feng et al., 2000). Nevertheless, the instant at which the leaflets start to close is approximately the same for both geometries.

Furthermore, for the asymmetrical geometry, the two leaflets show differences in movement. It can be understood that this asynchrony is triggered by the presence of the asymmetric geometry downstream of the valve (De Tullio et al., 2009; Hong et al., 2009). In the symmetrical geometries, there are no differences in movement between the two leaflets. The velocity magnitudes at different time levels for the three geometries are shown in Fig. 5.

Downstream of the valve, the three jet flows are clearly visible.

#### **5. Design challenges**

The ideal heart valve should have, among other things, a small drop in flow potential energy (i.e. small pressure drop over the valve and a high effective orifice area), a small retrograde flow, good hemodynamic properties, and high durability and safety in use. In the remainder of this section, each of these design challenges will be discussed in detail.

#### **Small regurgitation volume**

Mechanical valves are characterized by a significant amount of regurgitant flow. This retrograde flow is the sum of the closing volume and the leakage flow (Yoganathan et al., 2004). The closing volume is the amount of reverse flow needed to let the leaflets rotate to the closed position. The leakage flow occurs during diastole, when the leaflets are in the closed position and blood flows back through the gaps of the valve due to the large negative pressure gradient over the valve.

The regurgitation volume *Vreg* (in ml) is typically measured by percentage (%*reg*) of forward stroke volume *SVfwd* (in ml), as described in Verdonck et al. (2002):

$$\% \text{reg} = \frac{V\_{\text{reg}}}{V\_{\text{reg}} + SV\_{\text{fuel}}} \tag{17}$$

No turbulence model is used, thus implying laminar flow. A no-slip boundary condition is applied at the walls. The valve is initially set in the closed position. The moment of inertia of

The results of the simulations are depicted in Fig. 4 and Fig. 5. The angular positions of the leaflets are presented in Fig. 4(*a*), relative to the fully opened position. Therefore, 0 and 1

Although the valve leaflets open completely in the first geometry with the straight tube, the results show that this maximum opening position is not reached in both geometries with the sinuses of Valsalva. Such incomplete opening of the ATS Open PivotTM Standard Heart Valve in a divergent geometry is explained by the greater sensitivity of the leaflet movement to the flow field compared with other BMHV designs, since the leaflets extend farther in the flow downstream of the orifice than is the case in other valve designs (Feng et al., 1999). Therefore, the valve does not open completely in the divergent transvalvular flow caused by the enlargement of the sinuses of Valsalva because the leaflets tend to align with the

This incomplete opening also explains the difference in the closing phase between the geometries, since the valve in the sinuses of Valsalva geometries is closed sooner. This is because the leaflets in the straight tube reach the completely open position and therefore need to rotate over a greater angle in order to close. Hence, they have a greater closing volume (Feng et al., 2000). Nevertheless, the instant at which the leaflets start to close is

Furthermore, for the asymmetrical geometry, the two leaflets show differences in movement. It can be understood that this asynchrony is triggered by the presence of the asymmetric geometry downstream of the valve (De Tullio et al., 2009; Hong et al., 2009). In the symmetrical geometries, there are no differences in movement between the two leaflets. The velocity magnitudes at different time levels for the three geometries are shown in Fig. 5.

The ideal heart valve should have, among other things, a small drop in flow potential energy (i.e. small pressure drop over the valve and a high effective orifice area), a small retrograde flow, good hemodynamic properties, and high durability and safety in use. In the remainder of this section, each of these design challenges will be discussed in detail.

Mechanical valves are characterized by a significant amount of regurgitant flow. This retrograde flow is the sum of the closing volume and the leakage flow (Yoganathan et al., 2004). The closing volume is the amount of reverse flow needed to let the leaflets rotate to the closed position. The leakage flow occurs during diastole, when the leaflets are in the closed position and blood flows back through the gaps of the valve due to the large negative

The regurgitation volume *Vreg* (in ml) is typically measured by percentage (%*reg*) of forward

*reg fwd*

(17)

*V reg V SV* <sup>=</sup> <sup>+</sup>

% *reg*

one rigid valve leaflet about its rotation axis is equal to 9.94E-9kg·m2.

streamlines. In the straight tube, however, the valve leaflets open completely.

refer, respectively, to the fully closed and fully open position.

approximately the same for both geometries.

**5. Design challenges** 

**Small regurgitation volume** 

pressure gradient over the valve.

Downstream of the valve, the three jet flows are clearly visible.

stroke volume *SVfwd* (in ml), as described in Verdonck et al. (2002):

Since the regurgitation lowers the net cardiac output and therefore enlarges the workload of the heart, its percentage should be minimized when designing heart valves.

Fig. 5. Velocity Magnitude Contours in m/s visualized on a longitudinal section perpendicular to the leaflet rotation axes, for the geometry with the straight tube (*left*), the symmetrically placed sinuses of Valsalva (*middle*) and the asymmetrical sinuses of Valsalva (*right*). The plots are taken at *t* = 0.025s (*a*), *t* = 0.125s (*b*), *t* = 0.250s (*c*), and *t* = 0.375s (*d*), represented by the vertical lines in Fig. 4

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 43

One of the most challenging aspects in valve design is the improvement of hemocompatibility. Because the cardiovascular system is a closed circulatory system, damage to the blood cells can accumulate each time blood passes through the valve, which leads to an increased risk for platelet activation and hemolysis of the red blood cells (erythrocytes). This could result

Blood cells can become damaged by non-physiological flow patterns or by contact with artificial materials (Paul et al., 2003). Therefore, to optimize blood flow through the valve, one should minimize the production of shear stresses and turbulence in the flow and avoid any flow stagnation and separation in the vicinity of the valve (Yoganathan et al., 2004). Moreover, the valve should be composed of (or coated with) biocompatible materials.

Structural wear and fatigue of the valve can deteriorate the valve material and can cause

Moreover, cavitation, which leads to erosion and pitting, can occur in BMHVs. Such cavitation-induced erosion and pitting were first related to BMHVs in the 1980s after observation of a series of severe leaflet escapes of the Edwards-Duromedics bileaflet valve (Mastroroberto et al., 2000; Johansen, 2004). It was seen that the hard, but brittle, pyrolytic carbon can become subject to cavitation-induced fatigue leading to a transverse fracture of

Later, leaflet escape in the same (but revised) bileaflet valve type was observed in mitral

Therefore, for a durable and failure-safe design, the structural wear and fatigue potential should be minimized along with a minimization of the formation of cavitation bubbles.

Although current aortic BMHVs show a higher EOA and a higher PI compared with other artificial heart valve designs, such as the caged ball valve, the tilting disc valve, and the (non)stented bioprostheses (Yoganathan et al., 2004), some major design challenges concerning the hemodynamics and the durability need to be resolved. Moreover, BMHVs have a larger regurgitant volume than other artificial valves (Yoganathan et al., 2004). In the remainder of this section, the hemodynamics and the cavitation of the BMHV are

When blood flows through artificial devices, blood particulates can become damaged and initiate a cascade of events leading to coagulation and the formation of thrombi. Therefore, when a BMHV is implanted, the patient is required to take lifelong anticoagulation therapies (Paul et al., 2003; Bluestein et al., 2004; Dasi et al., 2009; Morbiducci et al., 2009). Blood trauma can occur via several mechanisms, depending on the type of damaged blood cell. In the past, the fragmentation and damaging of erythrocytes was experimentally quantified, since this leads to hemolysis (Sutera et al., 1975; Paul et al., 2003). In recent years, however, platelet activation is believed to be the major underlying formation mechanism for thromboembolic complications in the flow past mechanical heart valves (Bluestein, 2004;

the leaflet or a leaflet fracture near the pivot mechanism (Klepetko et al., 1989).

(Hemmer et al., 2000; Mert et al., 2003) and aortic position (Christiansen, 2001).

discussed in detail, with special interest to the flow near the hinges.

in clinical complications such as thromboembolisms and valve stenosis.

**Good hemodynamic properties** 

Therefore, mostly pyrolytic carbon is used.

**Durability and failure-safe design** 

**Design challenges and BMHVs** 

**5.1 Blood damage and BMHVs** 

Morbiducci et al., 2009).

severe valve failure.

#### **Low transvalvular pressure gradient (TPG)**

The transvalvular pressure gradient (TPG) is related to the drop in pressure that is needed for viscous blood to flow through the valve (Yoganathan et al., 2004). This forward pressure difference Δ*pDoppler* (in mmHg) can be derived from a simplified form of the Bernoulli equation and the mean (Doppler) velocity *Vmean* (in m/s) during forward flow (Verdonck et al., 2002):

$$
\Delta p\_{Doppler} = \mathbf{4} V\_{mean}^2 \tag{18}
$$

The TPG is a measure of valve efficiency (Yoganathan et al., 2004). The workload on the left ventricle is directly related to the magnitude of this pressure difference, because when the magnitude of the TPG increases, an increasing systolic pressure in the left ventricle is needed to drive flow forward into the aorta. Therefore, the TPG should be minimized when designing valve prostheses (Yoganathan et al., 2004).

#### **High effective orifice area (EOA)**

Because the TPG is heavily dependent on the flow rate, another less flow-rate-dependent parameter is more commonly used as a measure for the quality of the valve, namely the effective orifice area (EOA).

Several formulations of the EOA (in cm2) can be defined, of which the Gorlin equation is mostly used (Baumgartner et al., 1992; Yoganathan et al., 1984). It is based on the fundamental hydraulic law of flow through an orifice and couples the flow through the valve orifice to the Doppler pressure gradient over the valve (Verdonck et al., 2002):

$$EOA = \frac{Q\_{fuel,mem}}{51.6\sqrt{\Delta p\_{Doppler}}}\tag{19}$$

where *Qfwd,mean* is the mean forward flow rate (in cm3/s). The constant 51.6 accounts for gravity and unit conversions (Baumgartner et al., 1992; Verdonck et al., 2002; Yoganathan et al., 2007).

The EOA is a measure for valve quality because it assesses the severity of the stenosis formed by the presence of the valve and thus the degree to which the prosthesis obstructs the blood flow (Yoganathan et al., 2004). A large EOA corresponds to a smaller pressure drop and thus to a smaller energy loss (Yoganathan et al., 2004). Therefore, when designing a valve prosthesis, the resulting EOA should be maximized.

Moreover, the EOA can be used to compare the efficiency of various valve designs, because it is a better index of valve function than is TPG alone (Zoghbi et al., 2009). However, to get a normalized value, the EOA should be made independent of the valve size by dividing it by the valve sewing ring area *Asew*. This results in the performance index (PI), which provides a measure of the valve's resistance characteristics normalized to valve size (Yoganathan et al., 2004, 2007):

$$PI = \frac{EOA}{A\_{sw}}\tag{20}$$

Therefore, the parameter PI can be used to fairly compare the efficiency of various artificial valves.

#### **Good hemodynamic properties**

42 Aortic Valve

The transvalvular pressure gradient (TPG) is related to the drop in pressure that is needed for viscous blood to flow through the valve (Yoganathan et al., 2004). This forward pressure difference Δ*pDoppler* (in mmHg) can be derived from a simplified form of the Bernoulli equation and the mean (Doppler) velocity *Vmean* (in m/s) during forward flow (Verdonck et

The TPG is a measure of valve efficiency (Yoganathan et al., 2004). The workload on the left ventricle is directly related to the magnitude of this pressure difference, because when the magnitude of the TPG increases, an increasing systolic pressure in the left ventricle is needed to drive flow forward into the aorta. Therefore, the TPG should be minimized when

Because the TPG is heavily dependent on the flow rate, another less flow-rate-dependent parameter is more commonly used as a measure for the quality of the valve, namely the

Several formulations of the EOA (in cm2) can be defined, of which the Gorlin equation is mostly used (Baumgartner et al., 1992; Yoganathan et al., 1984). It is based on the fundamental hydraulic law of flow through an orifice and couples the flow through the

,

*fwd mean*

*Doppler*

valve orifice to the Doppler pressure gradient over the valve (Verdonck et al., 2002):

*<sup>Q</sup> EOA*

51.6

where *Qfwd,mean* is the mean forward flow rate (in cm3/s). The constant 51.6 accounts for gravity and unit conversions (Baumgartner et al., 1992; Verdonck et al., 2002; Yoganathan et

The EOA is a measure for valve quality because it assesses the severity of the stenosis formed by the presence of the valve and thus the degree to which the prosthesis obstructs the blood flow (Yoganathan et al., 2004). A large EOA corresponds to a smaller pressure drop and thus to a smaller energy loss (Yoganathan et al., 2004). Therefore, when designing

Moreover, the EOA can be used to compare the efficiency of various valve designs, because it is a better index of valve function than is TPG alone (Zoghbi et al., 2009). However, to get a normalized value, the EOA should be made independent of the valve size by dividing it by the valve sewing ring area *Asew*. This results in the performance index (PI), which provides a measure of the valve's resistance characteristics normalized to valve size

> *sew EOA PI*

Therefore, the parameter PI can be used to fairly compare the efficiency of various artificial

2

*Do* 4 *ppler mean* Δ = *p V* (18)

*<sup>p</sup>* <sup>=</sup> <sup>Δ</sup> (19)

*<sup>A</sup>* <sup>=</sup> (20)

**Low transvalvular pressure gradient (TPG)** 

designing valve prostheses (Yoganathan et al., 2004).

a valve prosthesis, the resulting EOA should be maximized.

**High effective orifice area (EOA)** 

effective orifice area (EOA).

(Yoganathan et al., 2004, 2007):

al., 2002):

al., 2007).

valves.

One of the most challenging aspects in valve design is the improvement of hemocompatibility. Because the cardiovascular system is a closed circulatory system, damage to the blood cells can accumulate each time blood passes through the valve, which leads to an increased risk for platelet activation and hemolysis of the red blood cells (erythrocytes). This could result in clinical complications such as thromboembolisms and valve stenosis.

Blood cells can become damaged by non-physiological flow patterns or by contact with artificial materials (Paul et al., 2003). Therefore, to optimize blood flow through the valve, one should minimize the production of shear stresses and turbulence in the flow and avoid any flow stagnation and separation in the vicinity of the valve (Yoganathan et al., 2004). Moreover, the valve should be composed of (or coated with) biocompatible materials. Therefore, mostly pyrolytic carbon is used.

#### **Durability and failure-safe design**

Structural wear and fatigue of the valve can deteriorate the valve material and can cause severe valve failure.

Moreover, cavitation, which leads to erosion and pitting, can occur in BMHVs. Such cavitation-induced erosion and pitting were first related to BMHVs in the 1980s after observation of a series of severe leaflet escapes of the Edwards-Duromedics bileaflet valve (Mastroroberto et al., 2000; Johansen, 2004). It was seen that the hard, but brittle, pyrolytic carbon can become subject to cavitation-induced fatigue leading to a transverse fracture of the leaflet or a leaflet fracture near the pivot mechanism (Klepetko et al., 1989).

Later, leaflet escape in the same (but revised) bileaflet valve type was observed in mitral (Hemmer et al., 2000; Mert et al., 2003) and aortic position (Christiansen, 2001).

Therefore, for a durable and failure-safe design, the structural wear and fatigue potential should be minimized along with a minimization of the formation of cavitation bubbles.

#### **Design challenges and BMHVs**

Although current aortic BMHVs show a higher EOA and a higher PI compared with other artificial heart valve designs, such as the caged ball valve, the tilting disc valve, and the (non)stented bioprostheses (Yoganathan et al., 2004), some major design challenges concerning the hemodynamics and the durability need to be resolved. Moreover, BMHVs have a larger regurgitant volume than other artificial valves (Yoganathan et al., 2004).

In the remainder of this section, the hemodynamics and the cavitation of the BMHV are discussed in detail, with special interest to the flow near the hinges.

#### **5.1 Blood damage and BMHVs**

When blood flows through artificial devices, blood particulates can become damaged and initiate a cascade of events leading to coagulation and the formation of thrombi. Therefore, when a BMHV is implanted, the patient is required to take lifelong anticoagulation therapies (Paul et al., 2003; Bluestein et al., 2004; Dasi et al., 2009; Morbiducci et al., 2009).

Blood trauma can occur via several mechanisms, depending on the type of damaged blood cell. In the past, the fragmentation and damaging of erythrocytes was experimentally quantified, since this leads to hemolysis (Sutera et al., 1975; Paul et al., 2003). In recent years, however, platelet activation is believed to be the major underlying formation mechanism for thromboembolic complications in the flow past mechanical heart valves (Bluestein, 2004; Morbiducci et al., 2009).

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 45

conditions for cavitation can be reached in these vortex cores (Avrahami et al., 2000; Johansen, 2004). Finally, the formation of cavitation bubbles can also be augmented by the sudden stop of the valve leaflets as they impact the casing, often referred to as "water

It is believed that in BMHVs none of these mechanisms alone but solely their combined interaction can initiate and augment cavitation. Moreover, valves with very stiff leaflets, closing at high velocity and decelerating rapidly (due to the impact at the blocking mechanism of the casing), are more prone to cause cavitation than valves with flexible

Although the formation of cavitation bubbles in the blood flow through a BMHV is an extremely rare phenomenon, it is undesirable because the shockwaves produced during the implosion of the cavitation bubbles can create high-velocity microjets. These microjets can deteriorate the valve structures as well as the nearby blood cells, thus leading to severe

Cavitation-induced material deterioration of a BMHV was first related to BMHVs in the 1980s after several leaflet escapes of the Edwards-Duromedics bileaflet valve, as described above (Mastroroberto et al., 2000; Johansen, 2004). It was observed that the impingement of high-velocity microjets can cause erosion and pitting when it impacts the structural surfaces. Lee et al. (2002) reported the appearance of cavitation-induced erosion pits in regions where squeeze flow occurred immediately before valve closure. Moreover, the study indicated that the number of pits was closely related to the magnitude of the pressure drop caused by the

Each year, replacing failing aortic heart valves with mechanical heart valves saves

However, modern BMHVs are still far from perfect and still face major design challenges. Most of these design challenges involve the hemodynamic properties of the valve and are thus directly related to the blood flow. Therefore, a thorough understanding of the blood

Current numerical simulation techniques can provide such valuable information and are considered crucial for gaining insights into the blood flow and assessing the performance of

The numerical simulations discussed in this chapter illustrate the incomplete opening of the aortic ATS Open PivotTM Standard Heart Valve in a diverging geometry. Moreover, asymmetrically placed sinuses of Valsalva induce an asymmetrical movement of the valve

Sebastiaan Annerel gratefully acknowledges a BOF-grant (Special Research Fund) from

Akutsu, T. & Higuchi, D. (2001). Flow analysis of the bileaflet mechanical prosthetic heart

valves using laser Doppler anemometer: effect of the valve designs and installed

leaflets that are gently decelerated (Zapanta et al., 1998; Johansen, 2004).

thromboembolic complications and valve failures (Johansen, 2004).

hammer cavitation"(Lee et al., 2002; Johansen, 2004).

water hammer phenomenon.

thousands of human lives.

future valve prototypes.

**7. Acknowledgment** 

**8. References** 

Ghent University Association.

leaflets.

flow is required for design optimization.

**6. Conclusion** 

The amount of blood damage depends primarily on the cumulative effect of the magnitude and the duration (the exposure time) of the applied force. Critical values of both factors can be exceeded by flow-dependent and non-flow-dependent causes (Paul et al., 2003).

The non-flow-dependent causes are the contact with artificial surfaces, which can be eliminated by using biocompatible materials. Therefore, pyrolytic carbon is commonly used for BMHVs because of its strength and high durability (Chandran et al., 2010) and good biocompatible properties (Johansen, 2004).

The flow-dependent causes are believed to originate from non-physiological flow patterns. However, despite many research efforts, the exact mechanisms underlying these flowinduced thromboembolic complications are still poorly understood (Bluestein, 2004). It is believed that the presence of elevated shear stresses in the flow is the most important platelet-damaging effect. Three non-physiological flow patterns can be distinguished.

Firstly, the squeeze flows observed as leakage flow when the leaflets are completely closed during diastole are of specific interest. Some leakage flow near the hinges is beneficial because it washes out the hinge regions and prevents flow stagnation and the development of thrombosis. However, the high-flow-velocity gradients can become too large, which causes elevated shear stresses and related platelet activation (Yoganathan et al., 2004).

Secondly, regions of flow separation, recirculation, and stagnation promote the deposition of damaged blood cells and increase the formation of thrombi (Yoganathan et al., 2004).

Finally, the shear layers surrounding the leaflets and the wake also expose the platelets to elevated shear stresses and lead them towards entrapment in the shed vortices of the wake (Bluestein et al., 2002).

Several studies used numerical methods to calculate the accumulated platelet activation when the blood flows through a BMHV. In the past, the valve leaflets were kept in a fixed position throughout the cardiac cycle due to the computational cost (Bluestein et al., 2002; Alemu et al., 2007; Dumont et al., 2005, 2007). More recently, however, Morbiducci et al. (2009) combined numerical FSI simulations with a numerical blood damage model. They reported that platelet activation is lower at early systole than at late systole and that the spanwise vorticity has a greater influence on the activation of platelets than does the streamwise vorticity.

A critical value for accumulated shear stress above which platelets are considered activated is, for example, given by the Hellums activation threshold criterion, i.e. 3.5Ns/m2 (Hellums et al., 1987; Dumont et al., 2007).

### **5.2 Cavitation and BMHVs**

Cavitation is the formation of voids or gas bubbles in liquids caused by a local reduction of the pressure to below that of vapor pressure. However, as soon as the surrounding pressure increases above vapor pressure, the formed bubbles will rapidly implode, which produces devastating shockwaves in the surrounding fluid (Johansen, 2004).

Several fluid mechanisms can lead to a pressure drop below vapor pressure in the blood flow during the closing phase of a BMHV. Firstly, squeeze flows are considered a contributing factor to the initation of cavitation. Such squeeze jets are formed at the very instant before leaflet closing, when the blood between the leaflets and valve casing is accelerated through the narrowing gap. This creates a high velocity jet flow with a large pressure gradient. The pressure can locally fall below vapor pressure, thus leading to cavitation (Bluestein et al., 1994; Johansen, 2004). Secondly, large vortices can be shed from the leaflet tips during the closing of the leaflets and during regurgitation. Towards the core of these vortices, the pressure decreases and the flow velocity increases. Therefore, the conditions for cavitation can be reached in these vortex cores (Avrahami et al., 2000; Johansen, 2004). Finally, the formation of cavitation bubbles can also be augmented by the sudden stop of the valve leaflets as they impact the casing, often referred to as "water hammer cavitation"(Lee et al., 2002; Johansen, 2004).

It is believed that in BMHVs none of these mechanisms alone but solely their combined interaction can initiate and augment cavitation. Moreover, valves with very stiff leaflets, closing at high velocity and decelerating rapidly (due to the impact at the blocking mechanism of the casing), are more prone to cause cavitation than valves with flexible leaflets that are gently decelerated (Zapanta et al., 1998; Johansen, 2004).

Although the formation of cavitation bubbles in the blood flow through a BMHV is an extremely rare phenomenon, it is undesirable because the shockwaves produced during the implosion of the cavitation bubbles can create high-velocity microjets. These microjets can deteriorate the valve structures as well as the nearby blood cells, thus leading to severe thromboembolic complications and valve failures (Johansen, 2004).

Cavitation-induced material deterioration of a BMHV was first related to BMHVs in the 1980s after several leaflet escapes of the Edwards-Duromedics bileaflet valve, as described above (Mastroroberto et al., 2000; Johansen, 2004). It was observed that the impingement of high-velocity microjets can cause erosion and pitting when it impacts the structural surfaces. Lee et al. (2002) reported the appearance of cavitation-induced erosion pits in regions where squeeze flow occurred immediately before valve closure. Moreover, the study indicated that the number of pits was closely related to the magnitude of the pressure drop caused by the water hammer phenomenon.

## **6. Conclusion**

44 Aortic Valve

The amount of blood damage depends primarily on the cumulative effect of the magnitude and the duration (the exposure time) of the applied force. Critical values of both factors can

The non-flow-dependent causes are the contact with artificial surfaces, which can be eliminated by using biocompatible materials. Therefore, pyrolytic carbon is commonly used for BMHVs because of its strength and high durability (Chandran et al., 2010) and good

The flow-dependent causes are believed to originate from non-physiological flow patterns. However, despite many research efforts, the exact mechanisms underlying these flowinduced thromboembolic complications are still poorly understood (Bluestein, 2004). It is believed that the presence of elevated shear stresses in the flow is the most important platelet-damaging effect. Three non-physiological flow patterns can be distinguished. Firstly, the squeeze flows observed as leakage flow when the leaflets are completely closed during diastole are of specific interest. Some leakage flow near the hinges is beneficial because it washes out the hinge regions and prevents flow stagnation and the development of thrombosis. However, the high-flow-velocity gradients can become too large, which causes elevated shear stresses and related platelet activation (Yoganathan et al., 2004). Secondly, regions of flow separation, recirculation, and stagnation promote the deposition of damaged blood cells and increase the formation of thrombi (Yoganathan et al., 2004). Finally, the shear layers surrounding the leaflets and the wake also expose the platelets to elevated shear stresses and lead them towards entrapment in the shed vortices of the wake

Several studies used numerical methods to calculate the accumulated platelet activation when the blood flows through a BMHV. In the past, the valve leaflets were kept in a fixed position throughout the cardiac cycle due to the computational cost (Bluestein et al., 2002; Alemu et al., 2007; Dumont et al., 2005, 2007). More recently, however, Morbiducci et al. (2009) combined numerical FSI simulations with a numerical blood damage model. They reported that platelet activation is lower at early systole than at late systole and that the spanwise vorticity has a greater influence on the activation of platelets than does the

A critical value for accumulated shear stress above which platelets are considered activated is, for example, given by the Hellums activation threshold criterion, i.e. 3.5Ns/m2 (Hellums

Cavitation is the formation of voids or gas bubbles in liquids caused by a local reduction of the pressure to below that of vapor pressure. However, as soon as the surrounding pressure increases above vapor pressure, the formed bubbles will rapidly implode, which produces

Several fluid mechanisms can lead to a pressure drop below vapor pressure in the blood flow during the closing phase of a BMHV. Firstly, squeeze flows are considered a contributing factor to the initation of cavitation. Such squeeze jets are formed at the very instant before leaflet closing, when the blood between the leaflets and valve casing is accelerated through the narrowing gap. This creates a high velocity jet flow with a large pressure gradient. The pressure can locally fall below vapor pressure, thus leading to cavitation (Bluestein et al., 1994; Johansen, 2004). Secondly, large vortices can be shed from the leaflet tips during the closing of the leaflets and during regurgitation. Towards the core of these vortices, the pressure decreases and the flow velocity increases. Therefore, the

devastating shockwaves in the surrounding fluid (Johansen, 2004).

be exceeded by flow-dependent and non-flow-dependent causes (Paul et al., 2003).

biocompatible properties (Johansen, 2004).

(Bluestein et al., 2002).

streamwise vorticity.

et al., 1987; Dumont et al., 2007).

**5.2 Cavitation and BMHVs** 

Each year, replacing failing aortic heart valves with mechanical heart valves saves thousands of human lives.

However, modern BMHVs are still far from perfect and still face major design challenges. Most of these design challenges involve the hemodynamic properties of the valve and are thus directly related to the blood flow. Therefore, a thorough understanding of the blood flow is required for design optimization.

Current numerical simulation techniques can provide such valuable information and are considered crucial for gaining insights into the blood flow and assessing the performance of future valve prototypes.

The numerical simulations discussed in this chapter illustrate the incomplete opening of the aortic ATS Open PivotTM Standard Heart Valve in a diverging geometry. Moreover, asymmetrically placed sinuses of Valsalva induce an asymmetrical movement of the valve leaflets.

## **7. Acknowledgment**

Sebastiaan Annerel gratefully acknowledges a BOF-grant (Special Research Fund) from Ghent University Association.

## **8. References**

Akutsu, T. & Higuchi, D. (2001). Flow analysis of the bileaflet mechanical prosthetic heart valves using laser Doppler anemometer: effect of the valve designs and installed

State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 47

Browne, P., Ramuzat, A., Saxena, R. & Yoganathan, A.P. (2000). Experimental investigation

Butany, J., Ahluwalia, M.S., Munroe, C., *et al.* (2003). Mechanical heart valve prostheses: identification and evaluation (erratum). *Cardiovasc Pathol*, Vol. 12, pp. 322-344. Carreau, P. (1972). Rheological equations from molecular network theories. *Trans Soc Rheol*,

Causin, P., Gerbeau, J.-F. & Nobile, F. (2005). Added-mass effect in the design of partitioned

Choi, C.R. & Kim, C.N. (2009). Numerical analysis on the hemodynamics and leaflet

Christiansen, S., Scheld, H.H. & Hammel, D. (2001). Leaflet escape in a Tekna valve in aortic

Cross, M. (1965). Rheology of non-Newtonian fluids: a new flow equation for pseudoplastic

Dahl, S.K., Vierendeels, J., Degroote, J., Annerel, S., Hellevik, L.R. & Skallerud, B. (2010). FSI-

Dasi, L.P., Ge, L., Simon, H.A., Sotiropoulos, F. & Yoganathan, A.P. (2007). Vorticity

Dasi, L.P., Simon, H.A., Sucosky, P. & Yoganathan, A.P. (2009). Fluid mechanics of artificial

Degroote, J., Bruggeman, P., Haelterman, R. & Vierendeels, J. (2008). Stability of a coupling

Degroote, J., Annerel, S. & Vierendeels, J. (2010). Stability analysis of Gauss-Seidel iterations

De Hart, J., Peters, G.W.M., Schreurs, P.J.G. & Baaijens, F.P.T. (2000). A two-dimensional

De Hart, J., Peters, G.W.M., Schreurs, P.J.G. & Baaijens, F.P.T. (2003). A three-dimensional

De Tullio, M.D., Cristallo, A., Balaras, E. & Verzicco, R. (2009). Direct numerical simulation

Diniz dos Santos, N.D., Gerbeau, J.-F. & Bourgat, J.-F. (2008). A partitioned fluid-structure

Donea, J., Huerta, A., Ponthot, J.-P. & Ferran A.R. (2004). Chapter 14 Arbitrary Lagrangian-

*Biomed Eng*, Vol. 28, No. 1, pp. 39-47.

method. *ASAIO J*, Vol. 55, No. 5, pp. 428-437.

*Method Biomec*, DOI: 10.1080/10255842.2010.517200.

heart valves. *Clin Exp Pharmacol Physiol*, Vol. 36, pp. 225-237.

position. *Ann Thorac Surg*, Vol. 71, p. 761.

systems. *J Colloid*, Vol. 20, pp. 417–437.

Vol. 16, No.1, pp. 99–127.

*Fluids*, Vol. 19, 067105.

Vol. 36, pp. 103–112.

Vol. 622, pp. 259-290.

Vol. 197, pp. 1750-1761.

470-84699-2.

2224-2234.

263-271.

pp. 4506-4527.

of the steady flow downstream of the St. Jude bileaflet heart valve: A comparison between laser Doppler velocimetry and particle image velocimetry techniques. *Ann* 

algorithms for fluid-structure problems. *Comput Methods Appl Mech Engrg*, Vol. 194,

dynamics in a bileaflet mechanical heart valve using a fluid-structure interaction

simulation of asymmetric mitral valve dynamics during diastolic filling. *Comput* 

dynamics of a bileaflet mechanical heart valve in an axisymmetric aorta. *Phys* 

technique for partitioned solvers in FSI applications. *Comput Struct*, Vol. 86, pp.

in a partitioned simulation of fluid-structure interaction. *Comput Struct*, Vol. 88, pp.

flui-structure interaction model of the aortic value. *J Biomech*, Vol. 33, pp. 1079-1088.

computational analysis of fluid–structure interaction in the aortic valve. *J Biomech*,

of the pulsatile flow through an aortic bileaflet mechanical heart valve. *J Fluid Mech*,

algorithm for elastic thin valves with contact. *Comput Methods Appl Mech Engrg*,

Eulerian Methods, In: *Encyclopedia of Computational Mechanics, Volume 1: Fundamentals*. E. Stein, R. de Borst, J.R. Hughes (Eds.), John Wiley & Sons, ISBN:0-

orientations to the flow inside the simulated left ventricle. *J Artif Organs*, Vol. 4, pp. 113-125.


Alemu, Y. & Bluestein, D. (2007). Flow-induced platelet activation and damage

Aljassim, O., Svensson, G., Houltz, E. & Bech-Hanssen, O. (2008). Doppler-catheter

Annerel, S., Degroote, J., Claessens, T. & Vierendeels, J. (2010). Evaluation of a new implicit

Aoyagi, S., Arinaga, K., Fukunaga, S., Tayama, E., Kosuga, T. & Akashi, H. (2006). Leaflet

Aslam, A.K., Aslam, A.F., Vasavada, B. & Khan, I.A. (2007). Prosthetic heart valves: types and echocardiographic evaluation. *Int J Cardiol*, Vol. 122, pp. 99-110. Astorino, M., Gerbeau, J.-F., Pantz, O. & Traoré, K.-F. (2009). Fluid-structure interaction and

Avrahami, I., Rosenfeld, M., Einav, S., Eichler, M. & Reul, H. (2000). Can vortices in the flow

Bang, J.S., Yoo, S.M. & Kim, C.N. (2006). Characteristics of pulsatile blood flow through the

vessels: velocity and pressure of blood flow. *ASAIO J*, Vol. 52, pp. 234-242. Baumgartner, H., Khan, S.S., DeRobertis, M., Czer, L.S. & Maurer, G. (1992). Doppler

Bech-Hanssen, O., Caidahl, K., Wallentin, I., Ask, P. & Wranne, B. (2001). Assessment of

Bluestein, D., Einav, S. & Hwang, N.H.C. (1994). A squeeze flow phenomenon at the closing

Bluestein, D., Li, Y.M. & Krukenkamp, I.B. (2002). Free emboli formation in the wake of bi-

Bluestein, D. (2004). Research approaches for studying flow-induced thromboembolic

Borazjani, I., Ge, L. & Sotiropoulos, F. (2008). Curvilinear immersed boundary method for

in vivo and in vitro study. *J Thorac Cardiovasc Surg*, Vol. 122, pp. 287-295. Bentall, H. & De Bono, A. (1968). A technique for complete replacement of the ascending

the aortic valve position. *Am J Cardiol*, Vol. 2008, pp. 1383-1389.

mm valves. *Ann Thorac Surg*, Vol. 82, pp. 853-857.

113-125.

No. 9, pp. 677-688.

*Biomec*, In press.

Vol. 198, pp. 3603-3612.

38, pp. 93-97.

pp. 2275-2283.

1378.

pp. 65-80.

aorta. *Thorax*, Vol. 23, pp. 338-339.

*Biomech*, Vol. 35, pp. 1533-1540.

Vol. 227, No. 16, pp. 7587-7620.

orientations to the flow inside the simulated left ventricle. *J Artif Organs*, Vol. 4, pp.

accumulation in a mechanical heart valve: numerical studies. *Artific Organs*, Vol. 31,

discrepancies in patients with bileaflet mechanical prostheses or bioprostheses in

coupling algorithm for the partitioned fluid-structure interaction simulation of bileaflet mechanical heart valves. *IOP Conf Ser: Mater Sci Eng*, Vol. 10, 012124. Annerel, S., Degroote, J., Claessens, T., *et al.* (2011). A fast strong coupling algorithm for the

partitioned fluid-structure interaction simulation of BMHVs. *Comput Method* 

movement of the ATS Valve in the aortic position: unique behavior observed in 19-

multi-body contact: application to aortic valves. *Comput Methods Appl Mech Engrg*,

across mechanical heart valves contribute to cavitation? *Med Biol Eng Comput*, Vol.

curved bileaflet mechanical heart valve installed in two different types of blood

assessment of prosthetic valve orifice area. An in vitro study. *Circulation*, Vol. 85,

effective orifice area of prosthetic aortic valves with Doppler echocardiography: An

of a bileaflet mechanical heart valve prosthesis. *J Biomech*, Vol. 27, No. 11, pp. 1369-

leaflet mechanical heart valves and the effects of implantation techniques. *J* 

complications in blood recirculating devices. *Expert Rev. Med Devices*, Vol. 1, No. 1,

simulating fluid structure interaction with complex 3D rigid bodies. *J Comput Phys*,


State-Of-The-Art Methods for the Numerical Simulation of Aortic BMHVs 49

Küttler, U. & Wall, W.A. (2008). Fixed-point fluid-structure interaction solvers with dynamic

Lee, H., Yamamoto, K., Kudo, N., Shimooka, T., Mitamura, Y. & Yuhta, T. (2002).

Le Tallec, P. & Mouro, J. (2001). Fluid structure interaction with large structural displacements. *Comput Methods Appl Mech Engrg*, Vol. 190, pp. 3039-3067. Makhijani, V.B., Yang, H.Q., Dionne, P.J. & Thubrikar, M.J. (1997). Three-dimensional

Masson, C. & Rieu, R. (1998). Time-frequency analysis of the noise produced by the closing of artificial heart valves: an in vitro study. *Med Eng Phys*, Vol. 20, pp. 418-431. Mastroroberto, P., Chello, M., Bevacqua, E., Cirillo, F. & Covino, E. (2000). Duromedics

Mert, K., Ozkara, A. & Hatemi, A. (2003). Leaflet escape in a revised Edwards-Duromedics

Mok, D.P., Wall, W.A. & Ramm, E. (2001). Accelerated iterative substructuring schemes for

Morbiducci, U., Ponzini, R., Nobili, M., *et al*. (2009). Blood damage safety of prosthetic heart

Morsi, Y.S., Yang, W.W., Wong, C.S. & Das, S. (2007). Transient fluid-structure coupling for

Nobili, M., Passoni, G. & Redaelli, A. (2007). Two fluid-structure approaches for 3D

Nobili, M., Morbiducci, U., Ponzini, R., *et al*. (2008). Numerical simulation of the dynamics

Paul, R., Apel, J., Klaus, S., Schügner, F., Schwindke, P. & Reul, H. (2003). Shear stress

Penrose, J.M.T. & Staples, C.J. (2002). Implicit fluid-structure coupling for simulation of

Peskin, C. (1972). Flow patterns around heart valves: a numerical method. *J Comput Phys*,

Redaelli, A., Bothorel, H., Votta, E., *et al*. (2004). 3-D Simulation of the St. Jude medical

Reul, H., Vahlbruch, A., Giersiepen, M., Schmitz-Rode, Th., Hirtz, V. & Effert, S. (1990). The

Sezai, A., Umeda, T., Hata, M., *et al*. (2009). A transesophageal echocardiographic and cine-

cardiovascular problems. *Int J Numer Meth Fl*, Vol. 40, pp. 467-478.

validation. *J Heart Valve Dis*, Vol. 13, pp. 804-813.

valve using closing velocities. *J Artif Organs*, Vol. 5, pp. 193–199.

Examination of cavitation-induced surface erosion pitting of a mechanical heart

coupled fluid–structure simulation of pericardial bioprosthetic aortic valve

original prosthesis: what do we really know about diagnosis and mechanism of

instationary fluid-structure interaction. In: *Computational Fluid and Solid Mechanics*,

valves. Shear-induced platelet activation and local flow dynamics: a fluid-structure

simulation of a trileaflet heart valve using weak coupling. *J Artif Organs*, Vol. 10,

simulation of St. Jude Medical bileaflet valve opening. *J Appl Biomater Biom*, Vol. 5,

of a bileaflet prosthetic heart valve using a fluid-structure interaction approaches. *J* 

related blood damage in laminar Couette flow. *Artif Organs*, Vol. 27, No. 6, pp. 517-

bileaflet valve opening process: fluid-structure interaction study and experimental

geometry of the aortic root in health, at valve disease and after valve replacement. *J* 

fluoroscopic evaluation of an ATS prosthetic valve opening. *Surg Today*, Vol. 39, pp.

relaxation. *Comput Mech*, Vol. 43, pp. 61-72.

function. *ASAIO J*, Vol. 43, No. 5, pp. M387–M392.

K. Bathe (Ed.), Elsevier, pp. 1325-1328.

*Biomech*, Vol. 41, pp. 2539-2550.

Vol. 10, pp. 252-271.

*Biomech*, Vol. 23, pp. 181-191.

pp. 96-103.

pp. 49-59.

529.

300-305.

approach. *J Biomech*, Vol. 42, pp. 1952-1960.

leaflet escape? *Can J Cardiol*, Vol. 16, No. 10, pp. 1269-1272.

mitral prosthesis. *J Heart Valve Dis*, Vol. 12, No. 4, pp. 513-515.


Dumont, K., Segers, P., Vandenberghe, S., Van Nooten, G. & Verdonck, R. (2002).

Dumont, K., Vierendeels, J., Segers, P., Van Nooten, G. & Verdonck, P. (2005). Predicting

Dumont, K., Vierendeels, J., Kaminsky, R., Van Nooten, G., Verdonck, P. & Bluestein, D.

Feng, Z., Umezu, M., Fujimoto, T., *et al*. (1999). Analysis of ATS leaflet behaviour by in vitro

Feng, Z., Umezu, M., Fujimoto, T., Tsukahara, T., Nurishi, M. & Kawaguchi, D. (2000). In

Feng, Z., Nakamura, T., Fujimoto, T. & Umezu, M. (2001). Influence of valve size on the hydrodynamic performance of the ATS valve. *J Artif Organs*, Vol. 4, pp. 303-307. Grigioni, M., Daniele, C., D'Avenio, G., *et al*. (2004). Innovative technologies for the

artificial heart valve testing. *Expert Rev Med Devices*, Vol. 1, No. 1, pp. 81-93. Guivier, C., Deplano, V. & Pibarot, P. (2007). New insights into the assessment of the

Guivier-Curien, C., Deplano, V. & Bertrand, E. (2009). Validation of a numerical 3-D fluid-

Hellums, J., Peterson, D., Stathopoulos, N., Moake, J. & Giorgio, T. (1987). Studies on the mechanisms of shear-induced platelet activation, Springer-Verlag, New York. Hemmer, W.B., Doss, M., Hannekum, A. & Kapfer, X. (2000). Leaflet escape in a TEKNA and an original duromedics bileaflet valve. *Ann Thorac Surg*, Vol. 69, pp. 942-944. Hirt, C.W., Amsden, A.A. & Cook, J.L. (1997). An arbitrary Lagrangian-Eulerian computing

Hufnagel, C.A., Harvey, W.P., Rabil, P.J. & McDermott, T.F. (1954). Surgical correction of

Irons, B.M. & Tuck, R.C. (1969). A version of the Aitken accelerator for computer iteration.

Johansen, P. (2004). Mechanical heart valve cavitation. *Expert Rev Med Devices*, Vol. 1, No. 1,

Kaminsky, R., Kallweit, S., Weber, H.-J., Claessens, T., Jozwik, K. & Verdonck, P. (2007).

Kelly, S., Verdonck, P., Vierendeels, J., Riemslagh, K., Dick, E. & Van Nooten, G. (1999). A

Klepetko, W., Moritz, A., Mlczoch, J., Schurawitzki, H., Domanig, E. & Wolner, E. (1989).

Flow visualization through two types of aortic prosthetic heart valves using stereoscopic high-speed particle image velocimetry. *Artif Organs*, Vol. 31, No. 12,

three-dimensional analysis of flow in the pivot regions of an ATS bileaflet valve. *Int* 

Leaflet fracture in Edwards-Duromedics bileaflet valves. *J Thorac Cardiovasc Surg*,

structure interaction model. *J Biomech*, Vol. 40, pp. 2283-2290.

method for all flow speeds. *J Comput Phys*, Vol. 135, pp. 203-216.

doppler study. *Int J Artif Organs*, Vol. 25, pp. 783-790.

*Heart Valve Dis*, Vol. 14, pp. 393-399.

experiment. *J Artif Organs*, Vol. 2, pp. 46-52.

position. *J Artif Organs*, Vol. 24, No. 5, pp. 346-352.

measurements. *Med Eng Phys*, Vol. 31, pp. 986-993.

aortic insufficiency. *Surgery*, Vol. 35, pp. 673–683.

*Int J Numer Meth Eng*, Vol. 1, pp. 275-277.

*J Artif Organs*, Vol. 22, pp. 754-763.

Vol. 97, No. 1, pp. 90-94.

pp. 95-104.

pp. 869-879.

pp. 558-565.

Omnicarbon 21 mm aortic valve prosthesis: in vitro hydrodynamic and echo-

ATS Open Pivot Heart Valve performance with computational fluid dynamic. *J* 

(2007). Comparison of the hemodynamic and thrombogenic performance of two bileaflet mechanical heart valves using a CFD/FSI model. *J Biomech Eng*, Vol. 129,

vitro hydrodynamic characteristics among three bileaflet valves in the mitral

assessment of cardiovascular medical devices: state-of-the-art techniques for

prosthetic valve performance in the presence of subaortic stenosis through a fluid-

structure interaction model for a prosthetic valve based on experimental PIV


**Part 2** 

**General Consideration of Aortic Valve Disease** 

