Skip to main content

Advertisement

Multipoint 5D flow cardiovascular magnetic resonance - accelerated cardiac- and respiratory-motion resolved mapping of mean and turbulent velocities

Abstract

Background

Volumetric quantification of mean and fluctuating velocity components of transient and turbulent flows promises a comprehensive characterization of valvular and aortic flow characteristics. Data acquisition using standard navigator-gated 4D Flow cardiovascular magnetic resonance (CMR) is time-consuming and actual scan times depend on the breathing pattern of the subject, limiting the applicability of the method in a clinical setting.

We sought to develop a 5D Flow CMR framework which combines undersampled data acquisition including multipoint velocity encoding with low-rank image reconstruction to provide cardiac- and respiratory-motion resolved assessment of velocity maps and turbulent kinetic energy in fixed scan times.

Methods

Data acquisition and data-driven motion state detection was performed using an undersampled Cartesian tiny Golden angle approach. Locally low-rank (LLR) reconstruction was implemented to exploit correlations among heart phases and respiratory motion states. To ensure accurate quantification of mean and turbulent velocities, a multipoint encoding scheme with two velocity encodings per direction was incorporated. Velocity-vector fields and turbulent kinetic energy (TKE) were obtained using a Bayesian approach maximizing the posterior probability given the measured data. The scan time of 5D Flow CMR was set to 4 min.

5D Flow CMR with acceleration factors of 19 .0 ± 0.21 (mean ± std) and velocity encodings (VENC) of 0.5 m/s and 1.5 m/s per axis was compared to navigator-gated 2x SENSE accelerated 4D Flow CMR with VENC = 1.5 m/s in 9 subjects. Peak velocities and peak flow were compared and magnitude images, velocity and TKE maps were assessed.

Results

While net scan time of 5D Flow CMR was 4 min independent of individual breathing patterns, the scan times of the standard 4D Flow CMR protocol varied depending on the actual navigator gating efficiency and were 17.8 ± 3.9 min on average. Velocity vector fields derived from 5D Flow CMR in the end-expiratory state agreed well with data obtained from the navigated 4D protocol (normalized root-mean-square error 8.9 ± 2.1%). On average, peak velocities assessed with 5D Flow CMR were higher than for the 4D protocol (3.1 ± 4.4%).

Conclusions

Respiratory-motion resolved multipoint 5D Flow CMR allows mapping of mean and turbulent velocities in the aorta in 4 min.

Background

The hemodynamic assessment of blood flow through heart valves and in the great vessels using 4D Flow cardiovascular magnetic resonance (CMR) has received increasing attention and a number of important applications have emerged [1]. While standard parameters including volume flow and peak velocities have been widely used, a number of additional readouts have been suggested. Among these readouts, probing the intravoxel velocity standard deviation [2] holds considerable promise as it permits estimating turbulent contributions in blood flow [3,4,5].

A number of clinical studies have indicated the added value of measuring turbulent velocities alongside mean velocities using 4D Flow CMR [2, 6, 7], e.g. by the assessment of turbulent kinetic energy (TKE) which indicates how much energy is stored in turbulent flow and which has shown elevated values for patients with bicuspid aortic valves or a dilated aorta relative to healthy subjects [7]. Ongoing studies promise further insights into the clinical relevance of turbulent flow estimation [8].

A key challenge of encoding mean and turbulent velocities relates to achieving appropriate measurement sensitivity to both quantities. For the expected ranges of mean and turbulent velocities in-vivo, multiple velocity encodings per spatial direction are preferably obtained and combined during image reconstruction [4]. This approach, however, requires additional data and hence data undersampling is necessary to arrive at clinically acceptable scan times.

Significant scan time reductions in 4D Flow CMR have been achieved using combinations of parallel imaging and compressed sensing [9], k-t methods [10, 11], radial [12] or spiral [13] acquisitions. Various modifications to the way data are sampled and reconstructed have been described, e.g. [14,15,16,17,18,19]. In general, however, scans remain too long to be performed in breathholds. Moreover, data acquisition without respiratory motion compensation leads to decreased image quality [14]. Accordingly, navigator-based respiratory gating is typically employed [20]. Unfortunately, respiratory gating makes scan times unpredictable, as gating efficiencies vary among subjects and often even during a single scan. Clinically, this may lead to conflicts as scan slots have typically fixed durations.

Instead of using respiratory gating, respiratory-motion resolved imaging [21,22,23,24] has been proposed to address respiratory-motion related image artifacts. To this end, data are acquired continuously throughout the respiratory motion cycle and sorted afterwards during the image reconstruction task into discrete motion states based on a respiratory motion signal surrogate [23, 25]. Accordingly, data correlation not only along the cardiac phase dimension but also along the respiratory motion dimension can be exploited.

Data correlations lead to low-rank properties of the decoding problem to be solved in image reconstruction. Various studies have demonstrated that these low-rank characteristics are most efficiently exploited using local rather than global low-rankedness using patch-based decomposition of the multi-dimensional data [26, 27].

The objective of the present work was to implement and test respiratory-motion resolved Bayesian multipoint 5D Flow CMR for mapping both mean and turbulent velocities in the aorta with a fixed scan time.

Methods

Data acquisition and respiratory motion binning

The data acquisition and motion binning process is illustrated in Fig. 1. Data are sampled using a pseudo-radial Cartesian sampling pattern which maps radial spokes onto a Cartesian grid [28, 29]. Using the Golden angle principle [30, 31], spatial incoherence of temporally adjacent frames is ensured. To reduce eddy current effects and keep the acoustic noise level during the measurement to a minimum, the tiny Golden angle concept was employed to avoid large jumps in k-space [31].

Fig. 1
figure1

5D Flow data acquisition. a Data are sampled using a Cartesian pseudo-radial tiny Golden angle sampling pattern [33]. b Respiratory motion state detection is derived from profile d(kx, ky = 0, kz = 0), which is repeatedly acquired and processed using a combination of principal component analysis, lowpass filtering, and coil-clustering [36]. Each spoke is composed of 21 k-space profiles, one of them always through the k-space center . With Tr = 3.3. ms, respiratory motion is sampled at ca. 14 Hz. c Acquired data are sorted into four discrete respiratory motion states such that the acceleration factor for each motion state is similar

Data-driven respiratory motion detection [32] was performed as illustrated in Fig. 2a) using repetitive sampling of profiles through the k-space center d(kx, t) = d(kx, ky = 0, kz = 0, t) as part of the pseudo-radial Cartesian sampling pattern. Upon inverse Fourier transform of d(kx, t), the projection X(x, t) of the excited volume resolved along the frequency encode direction x and over time t was obtained, showing signal modulation depending on respiratory and cardiac motion. Given TR = 3.3 ms and with one in 21 readouts through the k-space center, a respiratory sampling rate of ca. 14 Hz was obtained. The signal of each coil channel was treated as a separate signal \( {\mathrm{X}}_{\mathrm{i}}\in {\mathbb{C}}^{N_x\times {N}_t};i\in \left\{1,2,\dots, nCoils\right\} \), with Nx being the number of samples in readout dimension and Nt the number of readouts through the k-space center. Only the magnitude was considered, as the main variation in signal phase was found to be related to the different velocity encodings which were alternated per heart beat (beat-interleaved) and which led to phase variations which were too large to be eliminated by lowpass filtering.

Fig. 2
figure2

Data driven respiratory motion detection. a Respiratory motion is extracted from the profiles through the k-space center d(kx, ky = 0, kz = 0). The temporal evolution of readouts through the k-space center shows a modulation with respiratory motion and cardiac motion. b The 2nd principal component provides the main variation of the signal. c A lowpass filter limits the spectral components to the expected frequency range of respiratory motion. d) Coil clustering is used to combine the signals from different channels and provides a final estimate of respiratory motion

Singular value decomposition was performed and the 2nd singular vector of abs(XTX) was selected, as it provided the predominant motion for each coil channel (Fig. 2b). To isolate respiratory motion, a lowpass filter was applied to limit the spectrum to frequencies below 0.6 Hz (Fig. 2c). Next, coil channel clustering [28] was performed to determine the predominant dynamics among all channels which were assumed to be the respiratory motion signal (Fig. 2d). Finally, the average over all these signals was calculated to provide the final respiratory motion signal surrogate (Fig. 2e).

Motion-resolved image reconstruction

Using the respiratory motion surrogate signal, the acquired data were distributed across discrete motion bins. A total of four motion bins was defined and the width of each motion bin was set to achieve a similar undersampling factor for each respiratory motion state and cardiac phase by setting the bin boundaries to the 0.25-, 0.5- and 0.75-quantile of the respiratory motion signal. For some combinations of respiratory motion state, heart phase and velocity encoding, this approach led to very low sampling rates. Therefore, a data sharing approach [33, 34] was used to fill up frames with acceleration factors higher than 20 with data from neighboring respiratory motion states.

Accordingly, the final data vector reads \( \boldsymbol{d}\in {\mathbb{C}}^{N_x{N}_c{N}_{hp}{N}_{k_v}{N}_{rs}} \), with Nx being the number of samples in the spatial frequency domain, Nc denoting the number of coils, Nhp heart phases, \( {N}_{k_v} \) velocity encodings and Nrs respiratory motion states. A locally-low-rank (LLR) [26, 27] approach was used to exploit correlations among heart phases (hp) and respiratory motion states (rs). Image-domain data for each velocity encoding \( {\boldsymbol{i}}_{kv}\in {\mathbb{C}}^{N_x{N}_{hp}{N}_{rs}} \) were reconstructed separately as:

$$ {\hat{\boldsymbol{i}}}_{k_v}=\underset{{\boldsymbol{i}}_{k_v}}{\ \mathit{\arg}\min }{\left\Vert \varOmega \mathcal{FS}{\boldsymbol{i}}_{k_v}-{\boldsymbol{d}}_{k_v}\right\Vert}_2^2+\lambda {\sum}_{b=1}^{N_b}{\left\Vert\ {\mathcal{R}}_b\left({\boldsymbol{i}}_{k_v}\right)\right\Vert}_{\ast } $$
(1)

with the undersampling operator Ω, Fourier transform \( \mathcal{F} \), coil sensitivities \( \mathcal{S} \), k-space data dkv and regularization weight λ. As illustrated in Fig. 3a), the operator \( {\mathcal{R}}_b \) selects the b-th out of Nb blocks of size nx × ny × nz in the image from all Nhp heart phases and Nrs respiratory motion states and transforms them into a Casorati matrix with dimensions nxnynz × NhpNrs. By penalizing the nuclear norm of this local Casorati matrix, low-rankedness is enforced locally.

Fig. 3
figure3

Image reconstruction using locally low rank approach followed by Bayesian multipoint unfolding. a) A locally low-rank approach is employed for each velocity encoding separately. The locally low-rank model divides the image into 3-dimensional patches. Each patch is reordered into local Casorati matrices for which a low rank is enforced by penalizing the nuclear norm. Compared to a global Casorati matrix, the values of the singular values decrease more rapidly. b) Following reconstruction of the individual velocity encodings, for each Cartesian direction the different velocity encodings kv are combined using a Bayesian multipoint approach. A Bayesian probability model [4] provides posterior probabilities for mean velocity v and intra-voxel standard deviation σ given the measured signal S. v and σ are chosen such that the posterior probability is maximized, providing maps of turbulent kinetic energy (TKE) and mean velocities

For comparison, respiratory-motion resolved data were also reconstructed by penalizing total variation (TV) along the respiratory motion states and heart phases [22]:

$$ {\hat{\boldsymbol{i}}}_{k_v}=\underset{{\boldsymbol{i}}_{k_v}}{\ \mathit{\arg}\min }{\left\Vert \varOmega \mathcal{FS}{\boldsymbol{i}}_{k_v}-{\boldsymbol{d}}_{\boldsymbol{kv}}\right\Vert}_2^2+{\lambda}_{hp}\left.T{V}_{hp}\Big({\boldsymbol{i}}_{k_v}\right)+{\lambda}_{rs}\left.T{V}_{rs}\Big({\boldsymbol{i}}_{k_v}\right) $$
(2)

where λhp and λrs denote the regularization weights along cardiac phases and respiratory motion states, respectively.

Multipoint velocity encoding and decoding

Assuming a Gaussian distribution of intra-voxel velocities, the following signal model of mean velocity \( \overline{\boldsymbol{v}} \) and fluctuating component with standard deviation σ [2] was used:

$$ S\left({\boldsymbol{k}}_{\boldsymbol{v}}\right)={S}_0{e}^{\frac{-{\boldsymbol{\sigma}}^2{k}_v^2}{2}}{e}^{-i{\boldsymbol{k}}_{\boldsymbol{v}}\overline{\boldsymbol{v}}} $$
(3)

where kv is related to the first gradient moment of a bipolar velocity encoding gradient by \( {\boldsymbol{k}}_{\boldsymbol{v}}=\gamma {\int}_0^Tt\boldsymbol{G}(t) dt \), with T being the time of application of gradient waveform G(t) and determines the encoding velocity as \( {v}_{enc}=\frac{\pi }{k_v} \).

The signal model implies a trade-off between a sufficiently high encoding velocity venc to avoid phase wraps [1], and a sufficiently high encoding efficiency of the fluctuating velocity components [5] such that venc~σ. It has been demonstrated that this trade-off may lead to insufficient sensitivity to fluctuating velocity components and hence a multipoint encoding scheme using two instead of one kv encoding point per axis was implemented [4]. Using the probability map \( P\left(\overline{\boldsymbol{v}},\boldsymbol{\sigma} |D,S\right) \), pixel-wise estimates of mean velocity \( \overline{\boldsymbol{v}} \) and intra-voxel standard deviation σ given the model S and the measured image data D were obtained [4]. Data acquired with different kv′ s were combined by choosing the values for \( \overline{\boldsymbol{v}} \) and σ maximizing the posterior probability of the measured data:

$$ \hat{\overline{\boldsymbol{v}}},\hat{\boldsymbol{\sigma}}=\arg \underset{\overline{v},\sigma }{\max }P\left(\overline{\boldsymbol{v}},\boldsymbol{\sigma} |D,S\right). $$
(4)

An illustration of multipoint acquisition and Bayesian decoding is provided in Fig. 3b).

In-vivo study

All in-vivo work was performed upon written informed consent of the subjects and according to local ethics regulations. Data of the ascending aorta of nine healthy subjects (30 ± 11 years, 6 male) were acquired on a 3 T CMR system (Ingenia; Philips Healthcare, Best, the Netherlands) using the proposed 5D Flow CMR and a standard, navigator-gated 4D Flow CMR approach [1] with a spatial resolution of 2.5 × 2.5 × 2.5mm3, TE = 3.3 ms, TR = 4.9 ms, flip angle = 8 and 25 cardiac phases. The different velocity encodings were alternated per heart beat. The image acquisition matrix was Nreadout = 100 ± 9, Nphase = 101 ± 8 and Nslice = 20 ± 1 (mean ± std). Using the proposed method, velocities were encoded with venc = 50 cm/s and 150 cm/s per axis and an additional venc = 0 cm/s measurement. Scan time was fixed to 4 min. After sorting the data into 4 respiratory motion states, acceleration factors were 19.0 ± 0.21 (mean ± std). On average the view-sharing approach led to 61% of the total number of acquired samples being shared among motion bins in the study cohort.

Standard 4D Flow CMR was recorded with a single venc (venc = 150 cm/s) per axis (4D Flow Ref). Two-fold acceleration using parallel imaging [35] and standard pencil-beam navigator gating on the diaphragm (5 mm gating window) was employed with the 4D Flow reference protocol. Accordingly, the effective scan time depended on the breathing pattern of the subjects. All measurements were performed with retrospective cardiac gating and a 28-channel receiver coil.

Prior to image reconstruction, the 28-channel data were compressed to 8 virtual channels [36]. The locally low-rank (LLR) reconstruction (Equation 1) was solved using an implementation of the fast iterative shrinkage-thresholding algorithm (FISTA) [37] provided with the Berkeley advanced reconstruction toolbox (BART) [38]. A patch size of nx = ny = nz = 8 was used for the locally low-rank constraint. The TV reconstruction (Equation 2) was performed with an adapted version of the code provided with [22]. The standard 4D Flow reference data were reconstructed with MRecon (GyroTools LLC, Zurich, Switzerland). Sensitivity maps were estimated from a separate scan using the ESPIRiT method [39].

The regularization hyperparameters λ in Equation 1 and λhp and λrs in Equation 2 were optimized for best agreement of velocity components in systole with data of the 4D Flow reference measurement:

$$ \hat{\boldsymbol{\lambda}}=\underset{\boldsymbol{\lambda}}{\mathit{\arg}\min }{\left\Vert {\boldsymbol{v}}_{recon}\left(\boldsymbol{\lambda} \right)-{\boldsymbol{v}}_{ref}\right\Vert}_1. $$
(5)

Equation 4 was minimized using the bayesopt function in MATLAB R2017b (The MathWorks, Natick, Massachusetts, USA) [40]. Resulting hyperparameters were λ = 0.01 for Equation 1 and λhp = 0.04 and λrs = 0.05 for Equation 2. For both optimizations, the number of iterations was set to 80.

Data analysis

Upon image reconstruction, concomitant gradient terms were corrected [41] and third-order background phase correction was applied [42, 43]. Segmentation of the aorta was performed using ITK-SNAP [44]. The same masks were used for 4D Flow reference data and data in expiration obtained with 5D Flow LLR and 5D Flow TV. Maximum intensity projections (MIP) of velocity magnitude and TKE maps in systole were calculated as:

$$ {\mathit{\operatorname{Im}}}_{MIP}\left(x,y\right)=\arg \underset{z}{\max}\left(\mathit{\operatorname{Im}}\left(x,y,z\right)\right). $$

In order to avoid noise at the vessel boundaries the segmentation masks were eroded by one pixel for MIP projections.

The normalized root-mean-square error (nRMSE) of the velocity magnitude in the segmented aorta was calculated according to:

$$ nRMSE=\sqrt{\frac{\sum_{ROI}{\left(\boldsymbol{v}-{\boldsymbol{v}}_{ref}\right)}^2}{\left| ROI\right|\ast \max \left({{\boldsymbol{v}}_{ref}}^2\right)}} $$

where |ROI| corresponds to the number of voxels in the segmentation mask.

TKE was calculated voxel-wise as:

$$ TKE=\frac{\rho }{2}\left({\sigma}_x^2+{\sigma}_y^2+{\sigma}_z^2\right) $$

where \( {\sigma}_j^2 \) denotes the variance of velocity deviations per direction j and ρ is the fluid density (here we assumed 1060 kg/m3 for blood).

Moreover, multi-planar reslicing along the aorta was conducted using VMTK (www.vmtk.org) and peak flow and peak through-plane velocities (i.e. the maximum velocity component in the corresponding cross-section of the aorta) were assessed for each plane. Peak through-plane velocities were determined upon application of a 3x3x3 median filter to reduce contributions from noise.

The effectiveness of respiratory-motion resolved 5D Flow CMR was assessed by comparing it to reconstructions using only data from the end-expiratory bin (EXP) and to reconstructions using no respiratory binning (NG).

To assess the performance of navigator gating in the standard 4D Flow CMR protocol, the positions measured by the navigator on the diaphragm were interpolated and each readout was attributed its corresponding navigator position. Then, mean and standard deviations of the acceptance window per heart phase were calculated. Due to the coarse temporal resolution of the navigator signal (one navigator per heart beat), linear interpolation was chosen in order to provide a lower bound on the width of the acceptance window for different cardiac phases.

Statistical analysis

Bland-Altman analysis [45] of peak flow and peak through-plane velocities was performed to assess the agreement of velocity fields obtained with 5D Flow LLR and the 4D Flow reference.

Results

Scan durations for 4D Flow CMR were 17.8 ± 3.7 min compared to 4 min for the 5D Flow protocol.

Figure 4 provides an analysis of the acceptance window for respiratory gating with one pencil beam navigator per cardiac cycle as being used in the 4D Flow reference measurements. If the navigator signal was within the acceptance window, data of all cardiac phases were accepted. This results in a wider acceptance window for later cardiac phases leading to corresponding uncertainty about the actual motion states and increasing image artifacts for later cardiac phases.

Fig. 4
figure4

Analysis of navigator-based respiratory gating in standard 4D Flow CMR. a) For an exemplary 4D Flow scan the diaphragm position measured by the pencil beam navigator over time is shown together with accepted data (red) and rejected data (yellow). When the navigator signal is within the acceptance window, data of all cardiac phases are accepted until the next navigator position is obtained. As a result, the effective acceptance window of the data varies as a function of heart phase and becomes wider towards diastole as shown in b)

Reconstructions exploiting respiratory motion states (5D Flow LLR) are compared to using only data in end-expiration (EXP LLR) and without gating (NG LLR) in Fig. 5. Magnitude images and TKE maps from 5D Flow LLR reconstructions show less artifacts compared to EXP LLR and NG LLR reconstructions.

Fig. 5
figure5

Comparison of 5D Flow LLR reconstructions in expiration relative to reconstructions using end-expiratory data only (EXP LLR) and without gating (NG LLR). 5D Flow LLR shows reduced aliasing artifacts in the magnitude, velocity and TKE maps. For both EXP LLR and NG LLR, residual motion artifacts are present in the magnitude images and increased noise is observed in the TKE maps

Magnitude and velocity magnitude images for the end-expiratory motion state obtained with 5D Flow LLR and 5D Flow TV are compared relative to the 4D Flow reference in Fig. 6. It can be seen that maximum intensity projections of velocity magnitudes are similar for all methods while TKE maps and magnitude images provided by 5D Flow TV show more residual aliasing when compared to 5D Flow LLR. The root mean square estimate (RMSE) of the reconstructed images was lower for 5D Flow LLR compared to 5D Flow TV for both magnitude and velocity magnitude (8.7 vs. 13.5% for and 7.9 vs 9.8%, respectively). Maps of TKE for all 9 subjects obtained with 5D Flow LLR and the 4D Flow reference are compared in Fig. 7. In general, TKE maps derived from the 4D Flow reference show a higher noise level compared to 5D Flow data.

Fig. 6
figure6

Results in systole comparing 5D Flow LLR with 5D Flow TV and standard 4D Flow. For 5D Flow TV and the 4D Flow reference residual aliasing is observed in the magnitude images. Moreover, the TKE maps obtained with 5D Flow TV and the 4D Flow reference show more noise compared to 5D Flow LLR. Exemplary slices in the ascending aorta (AA) and descending aorta (DA) show qualitatively similar results for through-plane velocities in the ascending aorta whereas in the descending aorta artifacts can be observed in the in-plane velocity components for 5D Flow TV

Fig. 7
figure7

Comparison of TKE maps in systole for the 4D Flow reference versus 5D Flow LLR. In the 4D Flow reference data increased noise compared to 5D Flow LLR is seen. Moreover, TKE maps from 4D Flow show high values in the descending aorta whereas TKE values are lower in the descending aorta for most cases with 5D Flow LLR

Exemplary slices in the ascending and descending aorta show through-plane velocity profiles to be in agreement among the 4D Flow reference and 5D Flow LLR whereas 5D Flow TV shows artifacts in the descending aorta.

Exemplary magnitude and velocity magnitude images for different respiratory bins reconstructed with 5D LLR are shown in Fig. 8 along with an exemplary respiratory curve. Bins in expiration are narrower as the duration of expiration was longer compared to inspiration. The displacement of the heart and aorta is clearly seen in both magnitude and velocity maps. While image quality is acceptable, there is degradation present from end expiration to end inspiration.

Fig. 8
figure8

a) Exemplary magnitude and velocity magnitude images for different respiratory bins reconstructed with 5D LLR. Image quality reduces from end expiration to end inspiration and artifacts can be observed in the descending aorta in the magnitude and velocity images reconstructed in end inspiration. b): The respiratory curve shows that the expiratory bins have a narrower range of motion than inspiratory bins

Figure 9 compares velocity fields, peak velocities and peak flow obtained for all 9 subjects. The nRMSE between the velocity magnitudes obtained with 5D Flow LLR and the 4D Flow reference was 8.9 ± 2.1%. On average, 5D Flow LLR shows good agreement of peak velocities and peak flow relative to the 4D Flow reference (peak velocities: 3.1 ± 4.4%, peak flow: −2.4 ± 6.9%).

Fig. 9
figure9

Quantitative analysis of results. Scan times, peak through-plane velocities and peak flow for 5D Flow LLR reconstructions are compared to the 4D Flow reference. Duration of the 5D Flow LLR acquisition is on average ca. 4.5 times lower than for 4D Ref. Bland-Altman analysis shows a mean difference of 3.1% and − 2.4% and limits of agreement of ±8.6% and ±13.6%

Discussion

In this work, a respiratory-motion resolved Bayesian multipoint 5D Flow CMR approach has been implemented based on pseudo-radial tiny Golden angle Cartesian sampling in conjunction with locally low-rank image reconstruction to map mean and turbulent velocities in the aorta in a fixed scan time of 4 min.

By exploiting data from all respiratory motion states, the duration of 5D Flow CMR becomes independent of the individual respiratory motion patterns of the subjects. In comparison to standard 4D Flow CMR protocols [1], 5D Flow CMR is about 4.5 times shorter on average. Whereas pencil-beam based navigator gating in standard 4D Flow CMR cannot ensure motion-free data across the entire cardiac cycle, as demonstrated in Fig. 4, repetitive data-driven respiratory motion detection and continuous data acquisition of 5D Flow CMR provide motion estimates throughout the cardiac cycle.

It has been demonstrated that the LLR model allows to improve reconstruction accuracy relative to using data from end expiration only, as shown in Fig. 5. Of note, low-rank models have been used to exploit correlations in dynamic CMR [46,47,48] and have already been proposed for reconstructing respiratory-motion resolved data. However, existing approaches confine the reconstruction to a subspace defined by basis functions along the respiratory motion dimension either by combining cardiac and respiratory motion in a Casorati matrix [21] or by means of low-rank tensor factorization [49, 50]. These approaches either require the acquisition of a separate scan to estimate a subspace or a subspace has to be estimated from undersampled k-space data, e.g. via low-rank tensor factorization [49]. The first option increases scan time and can lead to problems when there is a mismatch between both scans (e.g. bulk motion) while the second option adds complexity to the reconstruction as another hyperparameter needs to be optimized for subspace estimation prior to reconstruction of the actual image. In comparison, the LLR model proposed in the present work can be applied straight-forwardly by penalizing the nuclear norm of the data rearranged in a Casorati matrix.

When comparing the 5D Flow LLR reconstruction with the 5D Flow TV approach, both methods showed good agreement of velocity data relative to the 4D Flow reference. However, TKE maps and magnitude images exhibit less noise when using the 5D Flow LLR approach. The inferior performance of the TV approach is associated with the piecewise constant solutions favored by TV which in turn leads to smoothing of peak values. To reduce underestimation of peak velocities, hyperparameters may be determined for best agreement with reference velocities. However, this led to non-optimal results for magnitude images which therefore still showed a high degree of aliasing noise.

In-vivo, good agreement of peak velocities and peak flow values of 5D Flow LLR with data derived from standard 4D Flow CMR was found. A small bias towards higher velocities was observed for 5D Flow LLR compared to the 4D Flow reference. This can be related to previous findings that respiratory motion leads to blurring of the image [20, 51] which corresponds to spatial low-pass filtering and is therefore likely to reduce peak velocities.

A multipoint encoding scheme [4] was incorporated into the 5D Flow approach. Thereby the velocity vector field was encoded with different sensitivities to mean and fluctuating velocities to provide an accurate assessment of mean and turbulent velocities over a large dynamic range. The 7-point 5D flow LLR implementation yielded TKE maps as they can be expected for healthy, young volunteers, with moderate values of up to 300 J/m2 and turbulence mainly occurring in the region of the flow jet in the proximal aorta. In comparison, TKE maps derived from the 4D Flow reference with 4-point encoding appeared noisy, due to only using a single encoding velocity of 150 cm/s. Moreover, the 4D Flow reference showed TKE values in the descending aorta which were as high as in the ascending aorta. This is considered unrealistic as transient/turbulent flow typically occurs in the region downstream of the aortic valve and not in the descending aorta in healthy, young volunteers.

The acquisition time for 7-point encoded 5D Flow was fixed to 4 min. For 10-point encoding, as it has been used when assessing stenotic valves [7], this would correspond to an acquisition time of ca. 6 min for 5D Flow, compared to 17.2 ± 4.7 min for accelerated 10-point 4D Flow acquisitions [7].

Image quality of data in inspiration as shown in Fig. 8 was lower compared to data reconstructed in the end-expiratory motion state with 5D Flow LLR. This can be explained with inspiration taking a smaller fraction of the duration of the respiratory cycle than expiration. Moreover, the respiratory motion curve does not reach the maximum value in inspiration for every respiratory cycle, whereas it reaches an end-expiratory plateau in most cases. As respiratory motion states were defined to obtain similar acceleration factors in each motion state, this led to a wider range of motion in inspiratory states than in expiratory states.

The present work assessed mean and fluctuating velocity vector fields in healthy subjects only. When applying the proposed 5D Flow method in patients with aortic stenosis, an extension to 10-point encoding [7] needs to be implemented for sufficient dynamic range. Adding 3 more points would lead to an increase in scan time by ca. 42% and can thus still easily be integrated into clinical workflows.

While the present study was focused on turbulence encoding, the proposed technique can be readily applied for the acquisition of standard 4D flow CMR data in the aorta and in the heart. For different applications the scan time will change as a function of the size of the acquisition matrix, but respiratory motion is still expected to be accurately modelled as a low-rank problem. Therefore, fixing scan time independently of respiratory motion should still be feasible. Moreover, the suggested multipoint encoding approach provides improved accuracy for lower velocities, similar to other multi-VENC approaches [52, 53]. In addition, the protocol can be used to measure flow in different respiratory motion states to assess respiration dependent flow-patterns, e.g. the Fontan circulation [54, 55]. To obtain better results in inspiration, one might define bins with equal range of motion, thus leading to varying acceleration factors for different bins. In order to still obtain sufficient sampling rates in each motion state, one would therefore need to increase scan time beyond the current 4 min limit.

A limitation of the present study is the lack of a ground truth data. The 4D Flow reference data were acquired in a separate scan and consistency of flow conditions cannot be expected over the entire scan session. In addition, the 4D Flow reference data showed considerable respiratory-motion related artifacts in later cardiac phases since the pencil-beam navigator was played out right after R-wave detection and only once per cardiac cycle (acceptance rates varied between 42 and 72%). Accordingly, the quantitative comparisons with the 4D Flow data were limited to systole in order to avoid respiratory-motion corrupted cardiac phases as much as possible.

Finally, a drawback of both 5D Flow LLR and 5D Flow TV are considerably longer reconstruction times when compared to standard 4D Flow CMR (ca. 35 min for 5D LLR and ca. 60 min for 5D TV vs. ca. 3 min for conventional 4D Flow CMR on a workstation with two 14 Core Intel Xeon E5–2680 CPUs and 256 GB RAM). In the future, this relative disadvantage may be addressed by employing variational neural network based reconstructions and their deployment on dedicated hardware [56, 57].

Conclusion

Respiratory motion resolved multipoint 5D Flow CMR allows for breathing-pattern independent mapping of mean and turbulent velocities in 4 min. The reduction in scan time allows for integration of the sequence into standard clinical workflows. Further in-vivo studies are now warranted to assess the performance of the method in relevant patient populations.

Availability of data and materials

The data that support the findings of this study are available from the corresponding author upon reasonable request subject to restriction on use by the Ethics Committee of the Canton of Zurich. A demo script with exemplary data will be provided online upon acceptance of this manuscript.

Abbreviations

CMR:

Cardiovascular magnetic resonance

CS:

Compressed Sensing

EXP:

end-expiratory

FISTA:

fast iterative shrinkage thresholding

LLR:

Locally-low-rank

MIP:

Maximum Intensity Projection

MP:

Multipoint

NG:

non-gated

nRMSE:

normalized root-mean-square error

RMSE:

root mean square error

TKE:

Turbulent kinetic energy

TV:

Total Variation

VENC:

Velocity encoding

References

  1. 1.

    Dyverfeldt P, Bissell M, Barker AJ, Bolger AF, Carlhäll C-J, Ebbers T, et al. 4D flow cardiovascular magnetic resonance consensus statement. J Cardiovasc Magn Reson. 2015;17(1).

  2. 2.

    Dyverfeldt P, Sigfridsson A, Kvitting JPE, Ebbers T. Quantification of intravoxel velocity standard deviation and turbulence intensity by generalizing phase-contrast MRI. Magn Reson Med. 2006;56:850–8.

  3. 3.

    Dyverfeldt P, Gårdhagen R, Sigfridsson A, Karlsson M, Ebbers T. On MRI turbulence quantification. Magn Reson Imaging. 2009;27:913–22. https://doi.org/10.1016/j.mri.2009.05.004.

  4. 4.

    Binter C, Knobloch V, Manka R, Sigfridsson A, Kozerke S. Bayesian multipoint velocity encoding for concurrent flow and turbulence mapping. Magn Reson Med. 2013;69:1337–45.

  5. 5.

    Knobloch V, Binter C, Gülan U, Sigfridsson A, Holzner M, Lüthi B, et al. Mapping mean and fluctuating velocities by Bayesian multipoint MR velocity encoding-validation against 3D particle tracking velocimetry. Magn Reson Med. 2014;71:1405–15.

  6. 6.

    Dyverfeldt P, Hope MD, Tseng EE, Saloner D. Magnetic resonance measurement of turbulent kinetic energy for the estimation of irreversible pressure loss in aortic stenosis. JACC Cardiovasc Imaging. 2013;6:64–71.

  7. 7.

    Binter C, Gotschy A, Manka R, Kozerke S. Turbulent kinetic energy assessed by multipoint 4D flow MRI provides additional information relative to echocardiography for the determination of aortic stenosis severity. Circ Cardiovasc Imaging. 2017.

  8. 8.

    Ha H. Kvitting J-PE, Dyverfeldt P, Ebbers T. 4D flow MRI quantification of blood flow patterns, turbulence and pressure drop in normal and stenotic prosthetic heart valves. Magn Reson Imaging. 55 July 2018;2019:118–27. https://doi.org/10.1016/j.mri.2018.09.024.

  9. 9.

    Tariq U, Hsiao A, Alley M, Zhang T, Lustig M, Vasanawala SS. Venous and arterial flow quantification are equally accurate and precise with parallel imaging compressed sensing 4D phase contrast MRI. J Magn Reson Imaging. 2013;37:1419–26.

  10. 10.

    Giese D, Wong J, Greil GF, Buehrer M, Schaeffter T, Kozerke S. Towards highly accelerated Cartesian time-resolved 3D flow cardiovascular magnetic resonance in the clinical setting. J Cardiovasc Magn Reson. 2014;16:42.

  11. 11.

    Schnell S, Markl M, Entezari P, Mahadewia RJ, Semaan E, Stankovic Z, et al. K-t GRAPPA accelerated four-dimensional flow MRI in the aorta: effect on scan time, image quality, and quantification of flow and wall shear stress. Magn Reson Med. 2014;72:522–33.

  12. 12.

    Gu T, Korosec FR, Block WF, Fain SB, Turk Q, Lum D, et al. PC VIPR: a high-speed 3D phase-contrast method for flow quantification and high-resolution angiography. Am J Neuroradiol. 2005;26:743 LP–749 http://www.ajnr.org/content/26/4/743.abstract.

  13. 13.

    Sigfridsson A, Petersson S, Carlhäll C-J, Ebbers T. Four-dimensional flow MRI using spiral acquisition. Magn Reson Med. 2012;68:1065–73.

  14. 14.

    Bollache E, Barker AJ, Dolan RS, Carr JC, van Ooij P, Ahmadian R, et al. K-t accelerated aortic 4D flow MRI in under two minutes: feasibility and impact of resolution, k-space sampling patterns, and respiratory navigator gating on hemodynamic measurements. Magn Reson Med. 2018;79:195–207.

  15. 15.

    Valvano G, Martini N, Huber A, Santelli C, Binter C, Chiappino D, et al. Accelerating 4D flow MRI by exploiting low-rank matrix structure and hadamard sparsity. Magn Reson Med. 2017;78:1330–41. https://doi.org/10.1002/mrm.26508.

  16. 16.

    Kim D, Dyvorne HA, Otazo R, Feng L, Sodickson DK, Lee VS. Accelerated phase-contrast cine MRI using k-t SPARSE-SENSE. Magn Reson Med. 2012;67:1054–64. https://doi.org/10.1002/mrm.23088.

  17. 17.

    Santelli C, Loecher M, Busch J, Wieben O, Schaeffter T, Kozerke S. Accelerating 4D flow MRI by exploiting vector field divergence regularization. Magn Reson Med. 2016;75:115–25.

  18. 18.

    Rich A, Potter LC, Jin N, Ash J, Simonetti OP, Ahmad R. A Bayesian model for highly accelerated phase-contrast MRI. Magn Reson Med. 2016;76:689–701.

  19. 19.

    Kwak Y, Nam S, Akçakaya M, Basha TA, Goddu B, Manning WJ, et al. Accelerated aortic flow assessment with compressed sensing with and without use of the sparsity of the complex difference image. Magn Reson Med. 2013;70:851–8.

  20. 20.

    Dyverfeldt P, Ebbers T. Comparison of respiratory motion suppression techniques for 4D flow MRI. Magn Reson Med. 2017;78:1877–82.

  21. 21.

    Sigfridsson A, Wigström L, Kvitting JPE, Knutsson H. K-t2 BLAST: exploiting spatiotemporal structure in simultaneously cardiac and respiratory time-resolved volumetric imaging. Magn Reson Med. 2007;58:922–30.

  22. 22.

    Feng L, Axel L, Chandarana H, Block KT, Sodickson DK, Otazo R. XD-GRASP: Golden-angle radial MRI with reconstruction of extra motion-state dimensions using compressed sensing. Magn Reson Med. 2016;75:775–88.

  23. 23.

    Cheng JY, Zhang T, Alley MT, Uecker M, Lustig M, Pauly JM, et al. Comprehensive multi-dimensional MRI for the simultaneous assessment of cardiopulmonary anatomy and physiology. Sci Rep. 2017;7:6–9.

  24. 24.

    Feng L, Coppo S, Piccini D, Yerly J, Lim RP, Masci PG, et al. 5D whole-heart sparse MRI. Magn Reson Med. 2018;79:826–38. https://doi.org/10.1002/mrm.26745.

  25. 25.

    Bastkowski R, Weiss K, Maintz D, Giese D. Self-gated golden-angle spiral 4D flow MRI. Magn Reson Med. 2018;80:904–13.

  26. 26.

    Trzasko J, Manduca A, Borisch E. Local versus global low-rank promotion in dynamic MRI series reconstruction. In: Proceedings of the 19th annual meeting of ISMRM; 2011. p. 4371.

  27. 27.

    Zhang T, Pauly JM, Levesque IR. Accelerating parameter mapping with a locally low rank constraint. Magn Reson Med. 2015;73:655–61.

  28. 28.

    Zhang T, Cheng JY, Chen Y, Nishimura DG, Pauly JM, Vasanawala SS. Robust self-navigated body MRI using dense coil arrays. Magn Reson Med. 2016;76:197–205.

  29. 29.

    Cheng JY, Hanneman K, Zhang T, Alley MT, Lai P, Tamir JI, et al. Comprehensive motion-compensated highly accelerated 4D flow MRI with ferumoxytol enhancement for pediatric congenital heart disease. J Magn Reson Imaging. 2016;43:1355–68.

  30. 30.

    Winkelmann S, Schaeffter T, Koehler T, Eggers H, Doessel O. An optimal radial profile order based on the golden ratio for time-resolved MRI. IEEE Trans Med Imaging. 2007;26:68–76.

  31. 31.

    Wundrak S, Paul J, Ulrici J, Hell E, Rasche V. A small surrogate for the golden angle in time-resolved radial MRI based on generalized fibonacci sequences. IEEE Trans Med Imaging. 2015;34:1262–9.

  32. 32.

    Larson AC, Kellman P, Arai A, Hirsch GA, McVeigh E, Li D, et al. Preliminary investigation of respiratory self-gating for free-breathing segmented cine MRI. Magn Reson Med. 2005;53:159–68.

  33. 33.

    Rasche V, De Boer RW, Holz D, Proksa R. Continuous radial data acquisition for dynamic MRI. Magn Reson Med. 1995;34:754–61.

  34. 34.

    Zhang S, Block KT, Frahm J. Magnetic resonance imaging in real time: advances using radial FLASH. J Magn Reson Imaging. 2010;31:101–9.

  35. 35.

    Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42:952–62.

  36. 36.

    Buehrer M, Pruessmann KP, Boesiger P, Kozerke S. Array compression for MRI with large coil arrays. Magn Reson Med. 2007;57:1131–9.

  37. 37.

    Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm. Soc Ind Appl Math J Imaging Sci. 2009;2:183–202.

  38. 38.

    Tamir JI, Ong F, Cheng JY, Uecker M, Lustig M. Generalized magnetic resonance image reconstruction using the Berkeley advanced reconstruction toolbox. In: Proceedings of the ISMRM 2016 data sampling and image reconstruction workshop. 2016. p. 9660006. http://www.ismrm.org/workshops/Data16/.

  39. 39.

    Uecker M, Lai P, Murphy MJ, Virtue P, Elad M, Pauly JM, et al. ESPIRiT—an eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA. Magn Reson Med. 2014;71:990–1001.

  40. 40.

    Snoek J, Larochelle H, Adams RP. Practical Bayesian optimization of machine learning algorithms. Adv Neural Inf Process Syst 2012;:2951–2959. doi:1206.2944.

  41. 41.

    Bernstein MA, Zhou XJ, Polzin JA, King KF, Ganin A, Pelc NJ, et al. Concomitant gradient terms in phase contrast MR: analysis and correction. Magn Reson Med. 1998;39:300–8.

  42. 42.

    Busch J, Giese D, Kozerke S. Image-based background phase error correction in 4D flow MRI revisited. J Magn Reson Imaging. 2017:1–10. https://doi.org/10.1002/jmri.25668.

  43. 43.

    Walker PG, Cranney GB, Scheidegger MB, Waseleski G, Pohost GM, Yoganathan AP. Semiautomated method for noise reduction and background phase error correction in MR phase velocity data. J Magn Reson Imaging. 1993;3:521–30.

  44. 44.

    Yushkevich, Paul A., Yang Gao, and Guido Gerig. "ITK-SNAP: An interactive tool for semi-automatic segmentation of multi-modality biomedical images." 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, 2016.

  45. 45.

    Altman DG, Bland JM. Measurement in medicine: the analysis of method comparison studies. Stat. 1983:307–17.

  46. 46.

    Liang Z-P. Spatiotemporal imaging with partially separable functions. IEEE Int Symp Biomed Imaging. 2007;2:988–91.

  47. 47.

    Pedersen H, Kozerke S, Ringgaard S, Nehrke K, Kim W. K-t PCA: temporally constrained k-t BLAST reconstruction using principal component analysis. Magn Reson Med. 2009;62:706–16.

  48. 48.

    Otazo R, Candès E, Sodickson DK. Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components. Magn Reson Med. 2015;73:1125–36.

  49. 49.

    He J, Liu Q, Christodoulou A, Ma C, Lam F, Liang Z-P. Accelerated high-dimensional MR imaging with sparse sampling using low-rank tensors. IEEE Trans Med Imaging. 2016;0062(c):1–1. https://doi.org/10.1109/TMI.2016.2550204.

  50. 50.

    Christodoulou AG, Shaw JL, Nguyen C, Yang Q, Xie Y, Wang N, et al. Magnetic resonance multitasking for motion-resolved quantitative cardiovascular imaging. Nat Biomed Eng. 2018;2:215–26. https://doi.org/10.1038/s41551-018-0217-y.

  51. 51.

    Wang Y, Riederer SJ, Ehman RL. Respiratory motion of the heart: kinematics and the implications for the spatial resolution in coronary imaging. Magn Reson Med. 1995;33:713–9.

  52. 52.

    Lee AT, Bruce Pike G, Pelc NJ. Three-point phase-contrast velocity measurements with increased velocity-to-noise ratio. Magn Reson Med. 1995;33:122–6. https://doi.org/10.1002/mrm.1910330119.

  53. 53.

    Herment A, Mousseaux E, Jolivet O, DeCesare A, Frouin F, Todd-Pokropek A, et al. Improved estimation of velocity and flow rate using regularized three-point phase-contrast velocimetry. Magn Reson Med An Off J Int Soc Magn Reson Med. 2000;44:122–8.

  54. 54.

    Gewillig M. The Fontan circulation. Heart. 2005;91:839–46. https://doi.org/10.1136/hrt.2004.051789.

  55. 55.

    Jarvis K, Schnell S, Barker AJ, Rose M, Robinson JD, Rigsby CK, et al. Caval to pulmonary 3D flow distribution in patients with Fontan circulation and impact of potential 4D flow MRI error sources. Magn Reson Med. 2018; February 2017:1205–18.

  56. 56.

    Vishnevskiy V, Sanabria SJ, Goksel O. Image reconstruction via Variational network for real-time hand-held sound-speed imaging. In: International workshop on machine learning for medical image reconstruction; 2018. p. 120–8.

  57. 57.

    Hammernik K, Klatzer T, Kobler E, Recht MP, Sodickson DK, Pock T, et al. Learning a variational network for reconstruction of accelerated MRI data. Magn Reson Med. 2018;79:3055–71.

Download references

Acknowledgements

The authors thank Giuseppe Valvano for his valuable advice regarding low-rank modelling of 4D flow and Richard Droste for contributing to initial experiments.

Funding

Financial support by the Swiss National Science Foundation SNF grant 320030_149567 and by ETH PASC is gratefully acknowledged.

Author information

SK and JW conceived of the study. JW implemented the sequence and reconstructions. JW and HD performed data acquisition. All authors participated in revising the manuscript and read and approved the final manuscript.

Correspondence to Jonas Walheim.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Ethics Committee of the Canton of Zurich, Switzerland, and all subjects provided written informed consent. Imaging was performed at the Zurich University Hospital, Zurich, Switzerland. Anonymized data was analyzed at ETH Zurich with approval by the mentioned authority. Written, informed consent was obtained before the experiment according to ethics approval and institutional guidelines.

Consent for publication

Written informed consent was obtained from all subjects for publication of their individual details and accompanying images in this manuscript. The consent forms are held in the subjects’ notes and are available for review by the Editor-in-Chief.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • 4D flow MRI
  • Velocity mapping
  • Turbulent kinetic energy
  • Respiratory motion compensation
  • Cartesian Golden angle
  • Low-rank image reconstruction
  • Cardiovascular magnetic resonance