5D Flow Tensor MRI to Efficiently Map Reynolds Stresses of Aortic Blood Flow In-Vivo

Diseased heart valves perturb normal blood flow with a range of hemodynamic and pathologic consequences. In order to better stratify patients with heart valve disease, a comprehensive characterization of blood flow including turbulent contributions is desired. In this work we present a framework to efficiently quantify velocities and Reynolds stresses in the aorta in-vivo. Using a highly undersampled 5D Flow MRI acquisition scheme with locally low-rank image reconstruction, multipoint flow tensor encoding in short and predictable scan times becomes feasible (here, 10 minutes), enabling incorporation of the protocol into clinical workflows. Based on computer simulations, a 19-point 5D Flow Tensor MRI encoding approach is proposed. It is demonstrated that, for in-vivo resolution and signal-to-noise ratios, sufficient accuracy and precision of velocity and turbulent shear stress quantification is achievable. In-vivo proof of concept is demonstrated on patients with a bio-prosthetic heart valve and healthy controls. Results demonstrate that aortic turbulent shear stresses and turbulent kinetic energy are elevated in the patients compared to the healthy subjects. Based on these data, it is concluded that 5D Flow Tensor MRI holds promise to provide comprehensive flow assessment in patients with heart valve diseases.

Imaging is playing an increasing role in assessing the hemodynamic and structural consequences of aortic valve diseases 1 . Time-resolved volumetric mapping of blood flow velocities using 4D Flow MRI 2 offers insights into changes of mean and peak velocities 3 , flow displacement 4 , vorticity and helicity 5 , wall shear rates 6 and relative pressure gradients 7 . Besides the assessment of time-resolved velocity vector fields, the intensity of stochastic velocity fluctuations as encountered in transient and turbulent flows can be probed 8 .
In general, turbulence dissipates energy and increases resistance to flow, generating additional load for the cardiovascular system 9 . Moreover, the effective coefficient of friction in turbulent flows is higher compared to normal flow and hence shear forces acting on the formed elements in blood are accentuated, potentially leading to blood cell damage [10][11][12] .
It has been shown that by quantifying Turbulent Kinetic Energy (TKE), i.e. the energy stored in velocity fluctuations, important additional information is obtained relative to current clinical information in heart valve patients 13 . Beyond quantifying TKE, all components of the Reynolds stress tensor (RST) may be obtained using appropriate changes of the MRI pulse sequence design 14 , a concept that has been validated using simulation and simplified in-vitro experiments recently 15,16 . Such an approach may offer improved mapping of pressure gradients across heart valves and stenotic vessel sections [15][16][17] .
A key practical challenge to quantifying the RST in-vivo relates to the extended scan times required in order to encode velocity fluctuations along the minimum number of six non-collinear axes. In addition, the dynamic range of velocity fluctuations encountered in-vivo demands at least two measurements along each non-collinear axis 18 , leading to scan times well beyond clinically acceptable limits.
The objective of the present work was to develop an approach to efficiently map the RST and hence turbulent shear stresses in-vivo within clinically acceptable scan times. Our approach is based upon recent advances in compressed sensing and sparse recovery of respiratory-motion resolved 4D Flow MRI data, which we have presented previously 19 . Here we propose a framework to efficiently quantify velocities and the RST using a highly undersampled acquisition scheme with locally low-rank image reconstruction 20,21 and multipoint encoding per axis including Bayesian estimation of average velocity per voxel as well as intravoxel velocity standard deviations 18 . We term this approach 5D Flow Tensor MRI.
Using a total of 19 velocity encodings, 5D Flow Tensor MRI requires 10 min of scan time and hence enables data acquisition in a clinical setting. To demonstrate accuracy and precision of 5D Flow Tensor MRI, results of computer simulations based on previously collected in-vivo data and in-vitro particle tracking velocimetry of valvular flow are shown. In-vivo proof of concept of 5D Flow Tensor MRI is demonstrated on patients with a bio-prosthetic heart valve revealing elevated turbulent shear stresses and turbulent kinetic energy compared to healthy controls.

Results
MRI data acquisition and reconstruction. Figure 1 illustrates the 5D Flow Tensor MRI concept including data acquisition, multipoint encoding, data reconstruction and Bayesian processing. Data are sparsely sampled using a Cartesian golden angle trajectory and retrospectively sorted into discrete respiratory motion states and cardiac phases 19 . Each velocity encoding is reconstructed separately using a locally low-rank reconstruction approach. Velocities are encoded in six non-collinear directions using three velocity encodings per axis to cover the range of turbulence intensity and mean velocities for patients and healthy controls.

Distributions of and sensitivity to intravoxel standard deviations.
To make an appropriate choice of the number and strength of velocity encodings per spatial axis, the distribution of velocity intravoxel standard deviations (IVSD) 8 was compared based on retrospective 4D Flow MRI data of patients with moderate and severe aortic valve stenosis (N = 28) and healthy controls (N = 9) collected as part of a previous study 13 . As shown in Fig. 2a, ISVD reaches up to 0.8 m/s in patients, while peak ISVD values of 0.3 m/s are measured in healthy controls. Since the MR signal magnitude is non-linearly related to ISVD, velocity encodings per axis need to be distributed in a non-equidistant manner. As illustrated in Fig. 2c, a velocity encoding (VENC) of 0.5 m/s shows high sensitivity to IVSD in the healthy controls whereas a VENC of 1.50 m/s is optimal to probe IVSD in the aortic stenosis patients. Figure 2d illustrates the resulting uncertainty in IVSD quantification with noisy data. Using Monte-Carlo simulations, for each value of IVSD σ, 10 5 samples with additive white Gaussian noise were generated and mean and standard deviation of the IVSD estimates σ est were determined. In case σ is too high or too low, σ est decreases in accuracy. Moreover, values of σ for which the signal magnitude vanishes cannot be discerned and lead to a plateau in the plot. As can be seen, an encoding velocity of 0.5 m/s, which would cover the range of IVSD in healthy aortae, cannot discern elevated values in patients. To ensure an accurate estimate of IVSD over the entire observed range, a distributed encoding scheme with 0.5 m/s, 1.5 m/s and 4.5 m/s is proposed. The first two values cover the range of turbulence, whereas the latter value prevents aliasing in the mean velocity field. Spatial resolution and Signal-to-Noise requirements. The effect of different signal-to-noise ratios (SNR) and the impact of image resolution on TKE and maximum principal turbulent shear stress (MPTSS) quantification was assessed using data previously acquired with particle tracking velocimetry (PTV) 22 as summarized in Fig. 3. For low SNR, an increase in mean values is observed for MPTSS. For TKE, the average mean values remain stable for low values of SNR (1.7% increase at 20 dB) while an increase in standard deviation is observed for decreasing SNR (e.g. 6.8% increase at 20 dB). At an SNR of 30 dB, as estimated for the in-vivo scans, MPTSS is overestimated by 3.6% on average whereas TKE values show no relevant increase in mean value (0.2%). Figure 3b shows the impact of different image resolutions for an SNR of 30 dB. Exemplary images show an increase of MPTSS and TKE at the jet core for increased voxel sizes. For large voxel sizes, the distribution of MPTSS values is skewed towards higher values with a corresponding increase in mean values and standard deviation. At a resolution of 2.5 mm, as used for the in-vivo exams, MPTSS are overestimated by 15.9% on average. TKE distributions are also skewed towards higher values for large voxel sizes with an overestimation of 3.1% at 2.5 mm.
Accuracy and precision of TKE and MPTSS quantification were simulated in a Monte-Carlo simulation with 40 repetitions. Mean and standard deviation over the repetitions are provided in Table 1 for varying SNR at the highest resolution and Table 2 for different resolutions at an SNR of 30 dB respectively. Table 1 shows an increase in the random error for decreasing SNR. However, the random error on mean and standard deviation of the value distribution remains below 1% for all metrics. Table 2 shows the effect of increasing voxel sizes for a fixed SNR of 30 dB. For increasing voxel sizes, a systematic overestimation can be observed for all metrics. Moreover, mean and standard deviation of MPTSS and TKE distributions show a higher random error for increased voxel sizes. In-Vivo measurements. Flow in the aorta of two patients with a bioprosthetic aortic valve (65 yrs, female with a SJM Trifecta Aortic Valve TFGT-21A, 21 mm, and 80 years, female with an Edwards SAPIEN 3, 23 mm) and two healthy controls (26 yrs, female and 58 yrs, female) was acquired using the 5D Flow Tensor MRI approach on a clinical 1.5 T MRI system (Philips Healthcare, Best, The Netherlands) and a 5-channel receive array. Figure 4a shows exemplary results in a single slice for a patient and a healthy control (patient 65 yrs, female, and volunteer 26 yrs, female). The highest values of TKE and MPTSS can be seen downstream of the bio-prosthetic valve in the patient.

Discussion
This study has demonstrated in-vivo turbulent flow assessment using 5D Flow Tensor MRI in clinically feasible scan times for the first time. A multi-point encoding scheme was employed to probe the mean and fluctuating velocity components using non-collinear encoding directions, similar to concepts used in diffusion tensor imaging 23 . The approach permits, besides the assessment of time-resolved mean velocity vector fields, the quantification of Reynolds stresses and hence turbulent kinetic energy and turbulent shear stresses in-vivo.
Distributions of IVSD in the aortae of healthy subjects and patients with aortic valve disease were analyzed to choose velocity encodings (Fig. 2a,b). As illustrated in Fig. 2c,d, the choice of velocity encoding has considerable impact on the accuracy of IVSD quantification. This makes the choice of an appropriate encoding velocity crucial, when using a single encoding velocity per axis as in conventional 4D Flow MRI. To probe IVSD with an increased dynamic range, encoding of the RST was combined with a multipoint scheme 18 . In the present study, encoding velocities of 0.5 m/s, 1.5 m/s and 4.5 m/s were selected. As indicated in Fig. 2a, encoding velocities of 0.5 m/s and 1.5 m/s cover the expected range of IVSD. The additional velocity encoding at 4.5 m/s was used to avoid phase wraps in the reconstructed mean velocity fields. Of note, the particular choice of encoding velocities was made with respects to the range of observed IVSD, to prevent aliasing artifacts in the mean velocity fields, and to make echo times not too long. However, the encoding scheme yields further potential for optimization. In particular, the use of advanced phase unwrapping methods 24 might allow to leave out the highest VENC.
Simulation of the MRI acquisition and encoding process revealed an overestimation of TKE and MPTSS for large voxel sizes. The overestimation amounted to about 3.1% for TKE and to approximately 15.9% for MPTSS at the given acquisition resolution of 2.5 mm and at an estimated SNR of 30 dB. The impact of image resolution can www.nature.com/scientificreports www.nature.com/scientificreports/ be related to the assumption of Gaussian intra-voxel velocity distributions in the derivation of turbulence, which is not fulfilled for coarse image resolutions as shown in previous studies 25 .
The impact of SNR on quantification of TKE and MPTSS was found to be relatively low compared to the impact of resolution. Starting at low SNR values below 25 dB, an overestimation of MPTSS was observed whereas www.nature.com/scientificreports www.nature.com/scientificreports/ TKE estimates were robust even at lower SNR values. SNR was estimated at ca. 30 dB in this study. In this range, noise played only a minor role in the assessment of TKE and MPTSS.
As shown in the Monte-Carlo simulation, the error in real-world experimental conditions is mostly due to a loss in accuracy for reduced image resolutions, whereas the random fluctuations for repeated experiments is comparatively low. However, increasing image resolution would lead to a decrease in SNR and noise would start to compromise the assessment of turbulent quantities. Therefore, rather than increasing acquisition resolution, efforts to mitigate the effect of large voxel sizes by e.g. data assimilation approaches 26 are considered potential future options.
The feasibility of 5D Flow Tensor MRI to quantify distributions of MPTSS and TKE in patients with a bio-prosthetic valve relative to healthy controls has successfully been demonstrated. Distributions of TKE and MPTSS revealed distinct differences, while differences in mean velocity magnitudes were partly overlapping (Fig. 4b). In the healthy controls, the highest values of TKE and MPTSS were found near the vessel walls, which can be attributed to partial volume effects (there were also some differences between the two volunteers which can be related to the difference in age 27 ). In contrast, flow downstream of the prosthetic valves showed highest values of MPTSS and TKE in the proximal aorta, reaching values of up to 500 Pa and 600 J/m 3 , respectively. MPTSS values were found to be below the threshold of elevated risk of red blood cell damage which was estimated between ca. 600 Pa 10 and 800 Pa 28 . While mechanical heart valves have been associated with blood cell damage 12 , modern bio-prosthetic valves typically do not lead to complications 29 . An increase in shear stresses without reaching a critical level was therefore expected. It should, however, be noted that the implantation of bio-prosthetic valves is primarily indicated in the elderly population, while mechanical heart valves are preferred in younger patients. Accordingly, future work using 5D Flow Tensor MRI should include patients with mechanical heart valves to assess and compare TKE and MPTSS levels.
Of note, the fixed scan time of 10 minutes which was set for the in-vivo study was sufficient for all subjects examined in this study. However, in cases where patient geometry requires a much larger field of view, an increase in scan time might be required.
A limitation of the present study is that no ground truth data was available to assess the accuracy of the in-vivo scans. Accordingly, computer simulations were used to provide estimates of accuracy and precision. However, the simulations were based on PTV measurements with a resolution of 0.625 mm. Thus, the reference data were already subject to some discretization error and availability of higher resolution ground truth data might show an even higher overestimation of turbulence. Another practical drawback relates to the long data reconstruction times (ca. 1.5 h to 2 h on a workstation with two 14 Core Intel Xeon E5-2680 CPUs and 256 GB RAM) which implies that data evaluation can only be performed after the scan session. Currently ongoing work is addressing this inconvenience by using variational neural networks 30 which have already been shown to perform compressed sensing reconstruction of standard 4D Flow MRI data in less than a minute 31 .
In conclusion, 5D Flow Tensor MRI provides comprehensive quantification of turbulent flow in clinically feasible scan times. Its ability to assess elevated TKE and MPTSS in-vivo has successfully been demonstrated. Efficient in-vivo turbulence quantification will contribute also to methods aiming at quantifying irreversible pressure loss downstream of heart valves and stenotic sections.   32 : In the one-dimensional case, assuming a Gaussian intra-voxel velocity distribution (IVSD) of variance σ 2 , the MR signal S(k v ) reads 8 : The magnitude of the complex-valued MR signal can be written as 14 : Analogous to diffusion tensor imaging 23 , the RST can be determined by encoding along six non-collinear directions and solving a system of linear equations. For six measurements along six different velocity encodings , the following set of equations is obtained: Accordingly, the elements of the RST can be calculated voxel-wise using the pseudoinverse: In this study, matrix H was designed according to: To mitigate the effect of non-linear encoding of the ISVD, a multipoint approach 18 was used to probe the velocity field at different encoding strengths. Figure 1b  , the velocities in the Cartesian coordinate system x y z can be written as: vn vn cart cart 1 1 A solution to this overdetermined system of linear equations is provided by the pseudo-inverse: Value range of intravoxel standard deviations. Datasets previously obtained in 9 healthy volunteers and 28 patients with aortic valve stenosis 13 were retrospectively analyzed to determine the range of IVSD occurring in the ascending aorta ( Fig. 2b shows exemplary slices with the corresponding region of interest). The data were acquired and reconstructed with multipoint acquisition and Bayesian reconstruction 18 . Values of VENC were 4.50, 1.50, and 0.50 m/s for patients and 2.00, 0.67, and 0.40 m/s for the healthy control group.
The ascending aorta was segmented manually. To assess the distribution of the ISVD for the two groups, the relative probability p( ) σ of different values of IVSD in the segmented region was calculated for each subject and the mean and standard deviation of σ p( ) were determined over the patient cohort and the healthy control group respectively. Spatial resolution and Signal-to-Noise requirements. MRI acquisitions with varying SNR and image resolution were simulated based on flow through a 64% stenosis measured with particle tracking velocimetry (PTV). Details on acquisition and processing of the PTV data can be found in 25,34 . The dynamic and kinematic viscosity were 5.82 × 10 −3 Pa and 4.85 × 10 −6 m 2 /s, respectively, and the fluid density 1200 kg/m 3 . The velocity-to-noise ratio was determined to be larger than 10 3 . The PTV data were mapped onto a voxel size of . × . × . 0 625 0 625 0 625 mm 3 .
Based on the PTV data, the MRI signal was calculated according to Eq. 4. Encoding velocities were 0.5, 1.5, and 4.5 m/s. To limit the effect of artifacts in the numerical study a median filter of size 3 was applied to the components of the RST.
To assess acquisition with different voxel sizes, the signal was transformed to k-space and sampled using a window function with a bandwidth inversely proportional to the desired downsampling rate. Complex-valued white Gaussian noise of different strength was added to the data to obtain the desired SNR Signal SD(Noise) (10) which was calculated over all velocity encodings.
In-Vivo measurements. In-vivo assessment of the RST was performed in two patients with bio-prosthetic aortic valves and two healthy controls on a 1.5 T MR system (Philips Healthcare, Best, The Netherlands). The study was approved by the Ethics Committee of the Canton of Zurich, Switzerland, and all subjects provided written informed consent. Data were acquired using a cardiac-and respiratory-motion resolved Cartesian tiny golden angle acquisition scheme 35,36 including the necessary velocity encodings for RST measurements. Acquisition and reconstruction of the data is illustrated in Fig. 1. During image reconstruction, data were sorted into four discrete respiratory motion bins. View sharing 37,38 among respiratory motion states was used to ensure a minimum acceleration factor of 35 for each frame. Scan parameters were: voxel size of 2.5 mm × 2.5 mm × 2.5 mm, 25 cardiac phases, multipoint flow tensor encoding with VENCs of 0.5 m/s, 1.5 m/s, and 4.5 m/s, TE/TR = 3.9 ms/6.0 ms and scan duration of 10 minutes compared to 71 minutes for a fully sampled scan (which could increase by a factor of ca. 2 when using respiratory navigator gating).
Prior to reconstruction, noise pre-whitening was performed based on noise statistics from a separate scan acquired without radio-frequency excitation. Data for each velocity encoding strength and direction were reconstructed separately with BART 39 , enforcing a locally low-rank model 20,21 along cardiac phases and respiratory motion states 19 . The signal estimate S kv is thus obtained by iterative minimization of the cost term with the undersampling operator Ω, Fourier transform  , coil sensitivities  and k-space data d kv . The operator  b selects the b-th out of N b blocks of size of size n n n 22 22 22 x y z × × = × × in the image from all N hp heart phase s and N rs respiratory motion states and transforms them into a Casorati matrix with dimensions × n n n N N x y z h p rs . The reconstruction favors solutions for which this local Casorati matrix is low-ranked by penalizing its nuclear norm. The regularization weight λ was set to λ = .
0 005. Both, block size and regularization weight were tuned for best agreement of magnitude images of the healthy control with a fully sampled reference measurement.
Following image reconstruction, only data in the expiratory motion state were considered for further processing. SNR in the measured data was determined using the pseudo-replica method 40 with 40 repetitions averaged over the ascending aorta and over the velocity encodings. Of note, approximate linearity is assumed with locally low-rank reconstructions. Accordingly, using Gaussian distribution of noise, the pseudo-replica method was considered the best approximation for SNR assessment. For quantitative evaluations of in-vivo data the ascending aorta was manually segmented using ITK-SNAP 41 and for the simulated data, the flow jet in was masked.

Statistical analysis.
Value-distributions of TKE, MPTSS, and velocity magnitude were investigated using a Gaussian kernel density estimate 42,43 . Moreover, mean and standard deviations of the distributions were assessed.
Accuracy and precision of TKE and MPTSS quantification were assessed in a Monte-Carlo simulation with 40 repetitions and mean and standard deviation over the experiment repetitions were determined. Ethics approval. 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.

Data availability
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.