Respiratory motion corrected 4D flow using golden radial phase encoding

To minimize respiratory motion artifacts while achieving predictable scan times with 100% scan efficiency for thoracic 4D flow MRI.


| INTRODUCTION
The 4D flow MRI allows for the assessment of time-resolved, 3D and 3-directional velocities. 1 Thanks to its increasing availability and robustness, it is to date being applied in several clinical studies such as valvular and cardiovascular diseases. 2 Nevertheless, the use of 4D flow MRI in a routine clinical setting or larger clinical studies remains hampered by its long examination times. Although advanced image reconstruction approaches such as parallel imaging or spatiotemporal acceleration techniques [3][4][5][6] have been able to significantly shorten acquisition times, respiratory motion can still lead to unacceptable long scan time. Respiratory gating is the most commonly used motion suppression technique for Cartesian 4D flow MRI. [7][8][9] For regular breathing patterns, scan efficiencies of 40% to 50% are achieved, but for irregular breathing patterns, this can drop to below 20%, 10,11 leading to scan times that are more than 5 times longer than the nominal scan duration.
To overcome the challenge of unpredictably long scan times, several studies have suggested to obtain 4D flow MRI under free-breathing without motion compensation. [12][13][14] Although flow quantification in major vessels did not vary significantly between ungated and gated scans, several studies have demonstrated that respiratory gating is needed for the extraction of wall shear stress and achieves reduced bias, smaller variability, and better agreement to 2D reference scans in phantom and in vivo studies for flow quantification. 15,16 Other approaches to reduce motion artifacts in 4D flow MRI include non-Cartesian acquisitions such as radial or spiral schemes. [17][18][19][20][21] Although minimizing breathing artifacts, non-Cartesian frequency encoding often lacks robustness to systematic errors such as phase errors due to eddy currents or blurring due to lengthened readout times. 22,23 Here, we propose to use motion-corrected image reconstruction for 4D flow MRI to minimize motion artifacts and ensure accurate flow quantification while ensuring 100% scan efficiency and predictable scan times for all subjects. The 4D flow MRI data are acquired with a golden radial phase encod ing (GRPE) scheme, which provides a respiratory selfnavigator and allows for the estimation of motion information from the data itself. [24][25][26][27][28] It allows for a straight-forward planning, the acquisition of an isotropic, large field of view (FOV) covering the heart and surrounding vessels, and runs at a predictable scan time of 15 min. Data were re-binned into cardiac phases based on the electrocardiograph (ECG) signal, while breathing motion information was extracted from the central k-space lines and corrected for during reconstruction yielding a breathing motion compensated 4D flow dataset (GRPE-MOCO). Flow values were then quantitatively compared with a thin-slab reference Cartesian 4D flow MRI of the aorta running at a similar total scan time (CART-REF). Additionally, GRPE-MOCO was compared with a dataset reconstructed without respiratory motion correction (GRPE-UNCORR). The method was evaluated in 9 healthy volunteers.

| Data acquisition
Data were acquired on a 1.5T system (Ingenia, Philips, Best, The Netherlands) using a 28-channel anterior and posterior coil array. Nine volunteers without known cardiovascular disease (4 female; age, 26.9 ± 6.6 years) underwent a 15-min 4D flow GRPE scan and a 4D flow CART-REF scan. The study was approved by the local ethics committee and written informed consent was obtained from all participants before the examination.

| GRPE
Cartesian frequency encoding was acquired in feet-head (FH) direction ( Figure 1A). Phase encoding was carried out along a radial grid in the 2D phase encoding plane k y -k z . 28 Phase encoding was, therefore, not carried out on a Cartesian grid but along non-Cartesian radial lines in the k y -k z plane. All phase encoding points along one radial line were obtained before increasing the angle by a golden radial increment of 111.24°. The GRPE scan was planned to cover the entire upper abdomen (i.e., heart, aorta, and pulmonary vasculature) with the following acquisition parameters: repetition time (TR)/ echo time: 3.8/2.0 ms, FOV: 250 × 250-270 × 250-270 mm 3 , acquired voxel size: 2.5 × 2.5 × 2.5 mm 3 , imaging matrix: 100-108 × 100-108 × 100, flip angle: 7°, readout bandwidth: 84.6 kHz. The scan was not synchronized to the ECG and no respiratory compensation was applied prospectively. The ECG signal, however, was recorded during the scan using a vector cardiogram and used to retrospectively sort the data into 24 cardiac phases. Flow encoding was performed in each direction consecutively using a 4-point symmetric acquisition scheme. 29 A velocity encoding value of 150 cm/s was chosen. Acquisition time was fixed to 15:03 min. rate), was reconstructed. Breathing motion was compensated by placing a pencil-beam navigator on the liver-diaphragm interface in combination with an acceptance window of 6 mm. A parallel imaging (SENSE) factor of 2 was used, resulting in a nominal scan time excluding navigator efficiency of 6:14 ± 0:32 min (range, 5:30-7:03 min).

| Respiratory self-navigator
A respiratory self-navigator was obtained from the central kspace line of each GRPE line. 28,30 Data acquisition was carried out such that one full GRPE radial spoke was acquired for one flow encoding direction, before repeating the same GRPE radial spoke for the next flow encoding direction. Therefore, a central k-space line was obtained every 100*TR = 380 ms. The overall signal intensities of the central k-space lines vary between different flow encodings; therefore, a respiratory self-navigator was calculated for each flow direction separately.

| Motion estimation
Based on the respiratory self-navigator signal, k-space data of each flow encoding direction was re-binned separately into 8 respiratory phases (RPh 1 -RPh 8 ) ( Figure 1B). The rebinning was carried out independently of the cardiac phase information; therefore, RPh 1-8 are averages over the cardiac cycles. Motion states were then reconstructed using autocalibrated iterative SENSE with temporal and spatial total variation regularization. 26,31 To estimate nonrigid motion fields (M r ) between respiratory states, a spline-based image registration algorithm was used. 32 A bending energy penalty ensured smooth M r allowing a robust motion estimation in the presence of residual undersampling artifacts. A normalized mutual information metric was used to maximize the similarity between images of different respiratory motion states. End-expiration of the flow-compensated data set was used as reference motion state.

| Motion correction
Based on the recorded ECG signal, k-space data were retrospectively re-binned into 24 cardiac phases using a sliding window approach with a window overlap of 20% on each side ( Figure 1C). Retrospective arrhythmia rejection was applied to exclude data from heart cycles deviating more than 20% from the average cycle length. For temporal k-t regularization, low spatial and high temporal resolution data (training data) were extracted from the central 25% of the k-space. To minimize truncation artifacts, a 2D F I G U R E 1 Overview. A, For each flow encoding direction, 3D k-space data are obtained with the GRPE trajectory. The read-out direction FH (k x ) remains Cartesian, whereas both phase encoding directions (k y , k z ) are obtained along radial spokes with an angular increment of 111.24°. B, Based on a self-navigator signal obtained from the central k-space lines of each GRPE spoke, data are split into 8 respiratory motion states and motion fields (M) are obtained using nonrigid spline-based image registration (REG). C, M are then used in a motion-corrected k-t SENSE reconstruction to obtain the final motion-corrected 4D flow images GRPE-MOCO Gaussian filter was applied and training data were reconstructed using a non-Cartesian iterative SENSE approach. GRPE data were then reconstructed using a self-regularized autocalibrated iterative non-Cartesian k-t SENSE approach. 33,34 At each reconstruction iteration, every cardiac phase was respiratory-motion corrected using the motion fields M r yielding a respiratory motion-corrected 4D flow reconstruction (GRPE-MOCO). The motion corrected image reconstruction was based on the general matrix description proposed by Batchelor et al 35 where the encoding operator E describing the acquisition process is extended by a sum over all motion states N: C describes the coil sensitivity information for each coil c, F is the Cartesian Fourier Operator, G is the interpolation operator from Cartesian to non-Cartesian k-space, and S selects the k-space points acquired for each respiratory bin r. 27,28 During image reconstruction, an iterative conjugate gradient approach was then used to minimize the difference between the acquired k-space data K and the motion-corrected image I MOCO : k-t SENSE regularization was added to minimize undersampling artifacts.
Additionally, 4D flow images were reconstructed without respiratory motion correction or respiratory gating by combining all acquired data irrespectively of their corresponding breathing phase (GRPE-UNCORR). This reduces the encoding operator described in Equation 1 to A schematic overview of the data acquisition and binning into respiratory and cardiac motion states is given in Supporting Information Figure S1, which is available online.

| Data analysis
Respiratory motion amplitudes (RMA) were measured as the amplitude of the motion vectors between end-expiration and end-inspiration. RMA was assessed in three regions of interest: right hemi-diaphragm (i.e., where a respiratory navigator would be placed), root of the aorta, and aortic arch. The ratio between RMA of the right hemi-diaphragm and the root of the aorta was also calculated. On the GRPE-MOCO and GRPE-UNCORR data, SV were additionally evaluated based on the conservation of mass: SV in the main pulmonary artery (PA) was compared with the sum of SV in the right PA (RPA) and left PA (LPA). Similarly, the difference in SV in the AAo (Qs) and PA (Qp) was calculated. Based on the SV, the commonly used metric of Qp/Qs ratios were calculated.
A paired Wilcoxon test was used to evaluate statistical significance of the differences of the compared quantities for CART-REF versus GRPE-MOCO, CART-REF versus GRPE-UNCORR, and internal validation. A P-value < 0.05 was considered as statistically significant.

| RESULTS
Data were successfully acquired and reconstructed in all volunteers. An overall acquisition time of 14:46 ± 5:44 min and 15:03 min was achieved in the CART-REF and GRPE sequence, respectively. In the CART-REF sequence, the scan efficiency due to respiration and ECG arrhythmias was 47.7 ± 11.5% and a mean number of 21.7 ± 0.7 heart phases were reconstructed, resulting in a mean temporal resolution of 43.9 ± 7.1 ms. The GRPE sequence was reconstructed to a fixed number of 24 heart phases, resulting in a temporal resolution of 41.7 ± 6.6 ms. The overall undersampling factor for a single cardiac phase combining data from all breathing phases was 3-4 in the GRPE data. Supporting Information Figure S2 visualizes the distribution of k-space data in the 2D phase encoding plane k y -k z for respiratory bins and cardiac phases. Figure 2 shows coronal slices in end-expiration (RPh 1 ) and end-inspiration (RPh 8 ) of a respiratory resolved GRPE dataset used for motion estimation. The subtracted image depicts motion state differences in the heart and abdominal organs (Figure 2, red and blue arrows, respectively). After application of the nonrigid motion field (M 8 ), differences are clearly reduced. Animations of the respiratory resolved images for all volunteers are shown in Supporting Information Figure S3.

| Motion estimation
RMA obtained from the motion fields were 10 ± 5.3 (range, 4.4-21) mm at the right hemidiaphragm, 6.2 ± 3.6 (range, 3.1-14.3) mm at the root of the aorta and 3.7 ± 2.7 (range, 1.2-10.1) mm at the aortic arch. The ratio of RMA between right hemi-diaphragm and aortic root was 0.6 ± 0.1 (range, 0.5-0.7).   Figure 4C-F show Bland-Altman plots for the internal validation using SV differences between the PA and RPA + LPA and between PA (Qp) and AAo (Qs). Pulmonary volumes in the GRPE-UNCORR and GRPE-MOCO show a mean difference of 10.2 ± 16.6 mL and 1.5 ± 16.6 mL, respectively. Similarly, Qp and Qs in the GRPE-UNCORR data show a higher mean difference of 6.1 ± 10.3 mL (Qp/Qs = 1.08 ± 0.12; range, 0.94-1.37) as compared to the GRPE-MOCO data, with a mean difference of −0.6 ± 10.0 mL (Qp/Qs = 1.00 ± 0.12; range, 0.82-1.15). No significant differences were found in the internal validation. All compared values are summarized in Supporting Information Table S2. Figure 5 shows coronal magnitude images and streamlines in the pulmonary arteries for both reconstructions in 2 exemplary volunteers, one with stronger (left) and one with less respiratory-induced motion (right). Unlike volunteer 2, the magnitude images of volunteer 1 show higher edge sharpness and better visibility of small vessels in the GRPE-MOCO as compared to the GRPE-UNCORR data. In this particular example, streamlines in the pulmonary arteries also indicate that flow in smaller vessels show improved visualization in the GRPE-MOCO data of volunteer 1.

| DISCUSSION
The present work proposes a novel 4D flow acquisition technique which is straight-forward to plan, runs at a predictable scan time, covers the entire upper abdomen, corrects for breathing motion and uses 100% of the acquired data while running F I G U R E 2 Respiratory motion correction. End-expiration (RPh 1 ) and end-inspiration (Rph 8 ) image. Difference image (RPh 1 -Rph 8 ) shows changes of liver and spline (blue arrows) and heart anatomy (red arrows) due to breathing. Motion estimation provides estimation of nonrigid transformation between end-expiration and end-inspiration phase (M 8 ). After applying M 8 end-inspiration phase is correctly transformed to endexpiration motion state, which minimizes image differences at a similar scan time as a standard thin-slab aortic 4D flow sequence.
Flow quantification comparisons with standard aortic 4D flow acquisitions revealed differences in-line with previously described alterations attributable to physiological changes. 15,17,36 Figure 5 shows that respiratory motion can cause blurring of the individual flow encoded images, which then impairs the quantitative 4D flow estimation leading to flow underestimation. The proposed nonrigid motion correction approach clearly improves image quality as well as flow quantification results. Increased accuracy with respect to CART-REF and Qp/Qs ratio was observed for the GRPE-MOCO data. GRPE-UNCORR data show a higher range of values for Qp/Qs in healthy volunteers of up to 1.37 compared with a maximum value of 1.15 using GRPE-MOCO. A Qp/Qs ratio > 1.5 can be used for classification of significant shunts 37 ; therefore, GRPE-MOCO could reduce the risk of false classification. The analysis of PF shows significant underestimation in the GRPE-UNCORR data and nonsignificant differences in the GRPE-CORR data. The effect of noncompensated motion on PF is in-line with a previously study using computational fluid mechanics. 16 In contrast to a typically respiratory-gated Cartesian acquisition, the proposed self-gating approach does not require the acquisition of a respiratory navigator, which assures the homogenous coverage of the entire cardiac cycle while the signal remains in the steady-state. The temporal resolution of the respiratory self-navigator obtained in this study was 380 ms, which was sufficient for all subjects. Nevertheless, for patients with very rapid breathing, a higher temporal resolution might be required. Although ECG information is recorded and used for re-binning, heart cycle variations do not influence data acquisition and arrhythmic heart cycles can be treated retrospectively. The amplitudes of respiratory motion vary strongly between different volunteers. The ratio between motion of the right hemi-diaphragm and aortic root was on average 0.6 as reported in previous studies. 38 Nevertheless, the range of the ratio was between 0.5 and 0.7, which means that a respiratory navigator placed on the right hemi-diaphragm with a window of 6.0 mm would have reduced the motion in the aortic root for 1 subject to 3.0 mm but for another only to 4.2 mm. The respiratory motion estimation presented here does not suffer from such subject specific variations and ensures accurate local motion minimization. Although residual motion in each respiratory bin may still be present, it remains small compared with the spatial resolution of the acquisition. For example, for the aortic root, the average residual motion amplitude is below 1 mm. Additional intra-bin motion correction could, however, be applied to further reduce this remaining motion. 39 In combination with the proposed motion correction approach, all acquired GRPE data were used for the final 4D flow reconstruction. Therefore, the predictability of scan time is independent of irregularities or changes in the breathing pattern during the scan. As unpredictable scan times and scan aborts due to low acceptance efficiencies and breathing position drifts currently hamper the wider use of 4D flow, the herewith presented approach is expected to alleviate this hurdle and make 4D flow easier integrable in clinical routine protocols. Even in the current study in which only healthy subjects were scanned, scan times for the respiratory-gated Cartesian acquisition showed variations from 9 to 28 min, which is in agreement to previous studies that have shown navigator efficiencies below 20%. 10,11 For the iterative non-Cartesian k-t SENSE reconstruction, the number of iterations and the weight of the k-t regularization were optimized based on one dataset. Supporting Information Figure S4 shows flow curves for different reconstruction parameters in one subject. Based on this analysis, the number of iterations and the weights were set to 12 and 3, respectively. As scan parameters were left unchanged between volunteers, these parameters were left constant throughout all reconstructions. The retrospective separation of data leads to a random distribution of k-space points in each cardiac phase. The distribution of data in each cardiac phase is different between different flow encodings and volunteers but leads to similar undersampling artifacts. Due to the undersampling properties of GRPE, this iterative reconstruction technique could minimize these incoherent undersampling artifacts and ensure high quality 4D flow data. More advanced image reconstruction approaches such as compressed sensing 40 or k-t PCA 41 could be used to further improve image quality and allow for a reduction in scan time.
One potential limitation is the acquisition of a very large FOV. Although acquisition time is similar to a standard 4D flow acquisition of a sagittal slab covering the aorta, the added value of hemodynamic information in the entire region may be questioned. In this study, the acquisition of a square FOV with isotropic spatial resolution was prescribed due to implementation ease; however, a freely adaptable FOV could be implemented and beneficial for the assessment of specific regions. 42 In the current feasibility study, the straightforward planning and improved signal-to-noise ratio when using the large FOV was, however, preferred.
As compared to so-called pseudo-radial or pseudo-spiral Cartesian k-space sampling, 43-46 the used GRPE leads to a variable density of acquired phase encoding points located on non-Cartesian coordinates. This requires gridding operation and, therewith, longer reconstruction times.
The use of motion correction was shown to increase accuracy of flow quantities. However, the co-registration of images from different respiratory phases assumes invariant hemodynamics for different respiratory states. This assumption may not be valid in specific vessels such as the vena cavae or in specific anatomies such as Fontan patients who show a respiratory dependent hemodynamic. 47 In these cases, the current technique may be used to extract respiratorydependent 4D flow datasets (5D flow).

| CONCLUSIONS
The proposed GRPE 4D flow technique allows for motioncorrected 3D and 3-directional velocity quantification by exploiting 100% of the acquired data. A predictable scan time, straight-forward planning, and large geometrical coverage is achieved and may help integration into clinical applications.

CONFLICTS OF INTEREST
Dr. Kilian Weiss is an employee of Philips Healthcare. Dr. Daniel Giese is an employee of Siemens Healthcare since December 2018. All of Dr. Giese's input related to the submitted manuscript, however, was performed prior or outside of his duties at Siemens. All other authors have no conflict of interest to declare.