Deep variational network for rapid 4D flow MRI reconstruction

Phase-contrast magnetic resonance imaging (MRI) provides time-resolved quantification of blood flow dynamics that can aid clinical diagnosis. Long in vivo scan times due to repeated three-dimensional (3D) volume sampling over cardiac phases and breathing cycles necessitate accelerated imaging techniques that leverage data correlations. Standard compressed sensing reconstruction methods require tuning of hyperparameters and are computationally expensive, which diminishes the potential reduction of examination times. We propose an efficient model-based deep neural reconstruction network and evaluate its performance on clinical aortic flow data. The network is shown to reconstruct undersampled 4D flow MRI data in under a minute on standard consumer hardware. Remarkably, the relatively low amounts of tunable parameters allowed the network to be trained on images from 11 reference scans while generalizing well to retrospective and prospective undersampled data for various acceleration factors and anatomies.


Introduction
4D flow MRI provides spatio-temporally resolved quantification of blood flow and offers great potential for the assessment of cardiovascular disease, e.g. aortic valve stenosis, atherosclerosis, or vessel wall remodelling 1 . However, clinical adaptation of the method has been hampered by long exam times.
Many efforts have been dedicated to accelerate flow acquisition by exploiting redundancies in the data. Partial Fourier imaging 2 has been used for moderate acceleration 3 , but the underlying assumption of a slowly varying phase has been shown to be incorrect for 4D flow MRI 4 . Parallel imaging (PI) 5 , which exploits the spatially varying sensitivity of receiver elements in the coil array has become a standard for accelerated imaging, but undersampling rates are limited by noise amplification 6 . The advent of compressed sensing (CS) 7 has enabled acceleration of 4D flow MRI by acquiring only a subset of k-space data and exploiting prior information about data regularities during reconstruction 8,9,10,11,12,13 , with typical acceleration factors ranging from 5 10 up to 27 13 . In particular, the locally low rank (LLR) regularized reconstruction 14 has been a successful technique, which iteratively balances the data fidelity cost and the singular norm of a patch matrix stacked over cardiac phases (see Methods for details). However, iterative reconstruction methods as used in CS increase reconstruction times considerably, implying that evaluation of 4D flow MRI data will typically happen when the subject has already been moved out of the scanner.
In recent years, deep neural networks have gained increasing popularity in MR image reconstruction. In the training stage, the neural network learns abstract features from a set of scans. After training, newly acquired data are reconstructed with very little computational effort by inference with the learned weights. This reduction in reconstruction times can facilitate the use of accelerated imaging methods in clinical practice. Moreover, reconstruction results can be superior to traditional CS methods 15,16 . Some approaches discard concepts of iterative image reconstruction altogether, e.g. by learning end-to-end mappings from k-space to image space 17 . As a downside, such networks usually require abundant amounts of high quality training data which is not available for high-dimensional flow MRI. Model-based neural reconstruction networks can also be designed to replicate the behaviour of an iterative reconstruction by interlacing nonlinear convolutional filters with an operation that enforces closeness of the current image estimate to the acquired data 15,18,16,19 , similar to the data fidelity step in an iterative shrinkage-thresholding algorithm 20 . A recent study 21 showed that neural network architectures which incorporate such an operation generalize better to different undersampling rates. In contrast, generic architectures which are solely based on convolutional layers can even lead to deteriorated image quality when the undersampling factor is increased, although one would expect the reconstruction result to improve when more information is available. The adversarial approach for training MR reconstruction networks 22,23 is usually aimed at improving perceptual reconstruction quality, such as image sharpness. Typically, this is achieved at the expense of reconstruction nRMSE (normalized root-mean-square error) 24 , which is critical for flow quantification.
In this work, an approach based on the idea of deep variational neural networks 15 is implemented for rapid 4D flow reconstruction, wich is referred to as FlowVN hereafter.. For comparsion, the 3D variational network (VN) architecture as presented in 15 was adapted for 4D flow data by using 3D filter banks operating on subsets of xyzt dimensional data, yielding the model that we refer to as HamVN. The FlowVN twork architecture replicates 10 steps of an iterative image reconstruction, while allowing for learnable spatio-temporal filter kernels, activation functions, and regularization weights in each iteration. It is demonstrated that based on training performed with retrospectively undersampled data of healthy subjects, FlowVN can accurately reconstruct pathological flow in a stenotic aorta in 21 seconds. Moreover, an imaging study with healthy subjects demonstrates good agreement of reconstructions from prospective undersampling with reference measurement.   Figure 1: Breathing resolved 4D flow data acquisition. a, Data are sampled using a Cartesian pseudo-radial tiny Golden angle sampling pattern. b, Time-resolved data are binned to end expiration using a combination of principal component analysis, low-pass filtering, and coilclustering. c, Acquired data are sorted according to heart phase and velocity encoding. d, Datasets used during training and evaluation. To conduct prospective evaluation, 7 healthy volunteers underwent accelerated CS and reference PI acquisitions during the same scan session.

4D flow MRI reconstruction with FlowVN
The FlowVN architecture improves HamVN in the following ways: (i) linear activations are used instead of radial basis functions (RBF), (ii) the network is conditioned on the sampling rate, (iii) exponential weighting of intermediate layers is used as regularization, (iv) real and imaginary parts of the signal are filtered by shared weights, (v) momentum is considered during gradient descent (GD) unrolling, and, (vi) the data term allows tunable activation functions. The network is trained for a wide range of acceleration factors by allowing acceleration dependent weighting of data consistency and filtering steps. As illustrated in Fig. 1a, for each velocity encoding direction, the k-space data is acquired using a Cartesian Golden angle sampling strategy, yielding variable density undersampling patterns in k-space. The signal of a total of 28 physical coils is compressed into 5 virtual coils via clustering 25 . The samples are then sorted into respiratory bins and data in the end-expiratory bin is used for reconstruction.
A deep variational network can be seen as a differentiable sequence of an unrolled numerical optimization scheme. To enable learning, such sequence is then relaxed by allowing tunable filter weights and activation functions. As described in Methods, we unroll K=10 steps of a gradient descent with momentum governed by a scalar α (k) : At each k-th layer P (k) the current complex-valued spatio-temporal image estimate is represented, while S (k) maintains a running average of update steps. The update step G (k) consists of the data consistency and regularization terms (see Methods and Supplementary Algorithm 1 for details), that are weighted according to the sampling rate M = 1/R (R is the acceleration factor) via tunable activation functions ϕ ur (M ), respectively. The data consistency term modulates the k-space data residual via an activation function and maps them back to the image space via a conjugate imaging operator. The regularization term at each layer contains 3D filters grouped into 4 banks, where each bank performs convolutions in 3 dedicated dimensions, namely xyz, xyt, xzt and yzt, therefore avoiding costly 4D convolutions. To avoid overfitting, we assume shared filters and activation functions that operate on real and imaginary components of the image. Note that both data and regularization terms do not assume correlations between real and imaginary parts of the signal, as highlighted in Fig. 2b.
The image estimate P (K) (B, Θ) of the final layer can be then seen as a function of the kspace samples B and network parameters Θ. To tune the network parameters Θ we minimize the layer-wise exponentially weighted 1 image reconstruction loss: over the retrospectively undersampled training dataset T , where P * is the ground truth image. Layer weighting is controlled by parameter τ ≥0: when τ ≈ 0, the reconstruction error is penalized equally across layers, therefore gradients of network parameters have lower variance during stochastic optimization, yielding faster convergence. On the contrary, when τ → +∞, only reconstruction at the final layer P (K) is minimized, which improves fitting accuracy on the training data. It is worth mentioning that τ controls the trade-off between training reconstruction residual and network regularity. Similarly to Landweber iterations 26,15 and deep supervision 27 , such implicit regularization penalizes irregular representations at intermediate layers and favors networks that can provide fast reconstruction. We propose to initialize τ with zero and then gradually increase it according to the training schedule (see Methods).  Figure 2: FlowVN architecture and training. a, Structure of FlowVN and its training strategy which uses a reduced field-of-view data. b, Single unrolled iteration block consisting of convolutional filtering, data fidelity, undersampling modulation and gradient descent momentum terms. Tunable parameters are highlighted in red. c, FlowVN at inference time yields 4D image reconstruction. d, Target training image error and velocity magnitude error in the aorta evaluated using training data. e, Data and gradient term weighting functions shown for each of 10 layers. f, Exemplary slices of 3D xyz filters and their corresponding activation functions at layer 5.
To demonstrate the validity of our approach, we note that the extracted velocity magnitude error in the aorta decreases simultaneously with the target reconstruction error during training, as shown in Fig. 2d, therefore indicating that the target 1 image error is a valid training surrogate. It can be seen from Fig. 2e that the regularization term is suppressed for lower acceleration factors R (higher sampling rate M). A subset of learned FlowVN parameters Θ is shown in Fig. 2f illustrating that learned convolutions perform direction-dependent filtering.

Retrospective and Prospective Evaluation
Reconstructed image magnitudes (for a single velocity encoding component), estimated velocity magnitudes and their errors of a healthy volunteer data for acceleration factor R=14 are illustrated in Fig. 3 for retrospectively undersampled data. Compared to CS-LLR and HamVN, the proposed FlowVN provides better reconstruction accuracy in terms of image magnitude and velocities. Scatter plot and correlation analysis further suggest that the velocity magnitude image estimated via FlowVN is in better agreement with ground truth. As shown in Supplementary Table 2 these observations extend to other acceleration factors R (6-22) as tested on 7 healthy volunteers. Fig. 4 indicates that FlowVN can accurately reconstruct the jet at the inlet section of the aorta for a patient with a pathological aortic valve.
The prospective undersampling acquisition results are reported in Fig. 5a Table 1 shows that the proposed FlowVN is 30 times faster than CS-LLR reconstruction.

Discussion
Practical learning-based image reconstruction can be traced back to dictionary learning methods 28,29 , where prior information is learned from image patches and then used as a sparsity-inducing regularizer for iterative reconstruction. Such an approach yields orders of magnitude longer reconstruction times, compared to modern deep learning approaches. A straightforward application of deep artificial neural networks has been suggested for learning reconstruction as a regression from the k-space 17 or zero-filled reconstructions 30 directly into the image space. Although tempting, such an approach might be unjustified, because the k-space and zero-filling artefacts have global dependence on image intensities. The advent of effective automatic differentiation systems 31,32 revitalized the idea of unrolling 33 and relaxing numerical schemes that can solve the original reconstruction problem. Following this approach, a number of deep neural network architectures were proposed 18,34,15,35 , that disentangle image acquisition and image prior models. Unrolling gradient descent reconstruction with tunable filters and activation functions yields the HamVN architecture proposed by Hammernik et al. 15 . One advantage of VNs is that, compared to other deep architectures, they employ relatively limited number of free parameters to tune, therefore they are less susceptible to overfitting.
In this work we have further developed the VN architecture 15,36,37 to accommodate high performance undersampled 4D flow reconstruction with limited training. Namely, we avoid exponential model complexity growth by avoiding 4D convolutions and by using separable 3D convolutions  that are shared for real and imaginary parts of the image. Furthermore, in contrast to the original HamVN 15 , we train our FlowVN for a wide range of undersampling factors by allowing the regularization term to depend on them. As illustrated in Fig. 2e, regularization scaled by ϕ ur decreases as more samples are available, while the data term ϕ ud stays constant for most of the layers. Such conditioning allows network training on a larger variety of artefacts and as it is necessary in practice, since for a given fixed acquisition time, the precise value of the undersampling factor is not known a priori and depends on breathing and cardiac motion patterns. We hypothesize that the wide range of acceleration factors which were used simultaneously to train the FlowVN provided a diverse collection of aliasing artefacts and enabled robust learning on a remarkably limited training set of 11 subjects. The exponential weighting of the layer-wise reconstruction loss (3) further regularized FlowVN parameters by penalizing the nonlinear behaviour presented in HamVN reconstructions. Supplemental Fig. 6 and Table 3 provide quantification of reconstruction accuracy effects attributed to the modifications proposed with FlowVN. In particular, modifications to the network architecture result in a model that can better adapt to data and yield higher accuracy for retrospectively undersampled experiments, while the proposed exponential weighting of the training loss improves accuracy of the prospective evaluation, which indicates better generalization ability. It is worth noting, that FlowVN has only 1% more tunable parameters compared to HamVN (c.f. Table 1), while improving reconstruction nRMSE by 23% (averaged over acceleration factors as given in Supplemental Fig. 6). We note that 4D flow MRI greatly benefits from using coil information during reconstruction (see Supplementary Table 3). Accordingly, comparison with single-coil reconstruction networks 18 has limited benefit.  The proposed FlowVN is a learning-based approach for reconstructing undersampled 4D flow MRI data in under a minute. For fixed reconstruction accuracy, FlowVN enables higher acceleration factors (12% improvement compared to CS-LLR image nRMSE at R=16) and does not introduce significant bias of peak flow estimates. The proposed reconstruction is 30 times faster than stateof-the art CS-LLR and 4.2 times faster than HamVN, due to using linear activation functions rather than RBFs, which requires computation of pairwise distances between control knots and image intensities. It is worth noting that FlowVN demonstrates high generalization ability, being able to preserve patient pathologies that were not present in the training data.

Compressed Sensing 4D Flow Reconstruction
PC MRI encodes flow velocity v(r, t) ∈ R 3 at spatial location r during cardiac phase t (1 ≤ t ≤ N t ) according to the following equation: where v enc is the velocity corresponding to a phase of ±π, y i , i = 0, . . . , 3 are the encoded velocity vector components. The four-point velocity encoding matrix is given as Therefore, flow velocity v can be calculated from the phase difference of reconstructed PC images ρ i . Let ρ it ∈ C Nv be a discretized image on a N x ×N y ×N z =N v grid corresponding to a cardiac phase t and velocity encoding i. Assuming Cartesian sampling on a regular N 1 ×N 2 ×N 3 =N s grid, the Fourier transform F ∈ C Ns×Nv , and N c coil sensitivity maps W k def = diag(c k ) ∈ C Nv×Nv define the spatial encoding operator E ∈ C NsNc×Nv : Considering a single velocity-encoded image sequence, let P ∈ C Nv×Nt and B ∈ C NsNc×Nt be stacked column-vectors of signals ρ and zero-filled k-space samples respectively, while M ∈ {0, 1} NsNc×Nt defines the undersampling mask. Iterative image reconstruction methods seek for a maximum a posteriori (MAP) solution defined by the following optimization problem: where the regularization term R enforces prior assumptions about image regularities. Herein we consider the local low-rank (LLR) regularization 14 to leverage image correlations among cardiac phases: where T i ∈ {0, 1} p 3 ×Nv is the corresponding p×p×p patch extraction operator, yielding N ptch patches, and · * is the nuclear norm. For LLR regularization the optimization problem (7) is convex and can be efficiently solved using operator splitting techniques such as the fast iterative shrinkage-thresholding algorithm (FISTA) 20 .

FlowVN Training
We employ a K=10 layer VN and perform 5·10 4 iterations of the ADAM algorithm (learning rate 10 −3 , β 1 =0.85, β 2 =0.98, batch size of 3) for training, during which we continually adjust τ = i opt ·10 −3 with i opt being the iteration number. On every layer, each 3D filter bank contains N f =8 filters of size n c =5 voxels. Activation functions ϕ{.} are parametrized by N knts =91 control knots with spacing ω=0.17: with gradients provided by the following formulas: The acquired zero-filled k-space B with undersampling mask M was normalized by M 1 B F . To enable backpropagation to be carried out with limited GPU memory, we employ spatiotemporal equivariance of the convolution and exploit the fact that k-space is fully sampled in the readout dimension k x for Cartesian acquisitions. Therefore, to draw a training sample, we perform random cropping of width w x and w t in dimensions x and t respectively and simulate Fourier encoding in k y k z dimensions as illustrated in Fig. 2(a). The network was implemented using the Tensorflow framework 32 . Fully-sampled and partial Fourier acquisition data from 11 healthy volunteers was used during training.

In Vivo Data Acquisition
As illustrated in Fig. 1, we used 11 subject for network training, and 7 healthy subjects and 1 patient for evaluation. All in-vivo work was performed upon written informed consent of the subjects and according to local ethics regulations.
Training datasets comprised 4D flow data measured in the aorta of 11 healthy subjects, 9 of them fully sampled and 2 acquired with partial Fourier 38 (factor 0.75×0.75).
For evaluation, data in the ascending aorta of 7 healthy subjects were acquired on a 3T Philips Ingenia system (Philips Healthcare, Best, the Netherlands) using a Cartesian 4-point referenced phase-contrast gradient-echo sequence with an encoding velocity v enc = 150 cm/s, a spatial resolution of 2.5×2.5×2.5 mm 3 , T E = 3.3 ms, T R = 4.9 ms, 25 cardiac phases and flip angle = 8 • . Exams for each of the 7 healthy subjects comprised a standard navigator-gated 2-fold accelerated parallel imaging 5 exam for reference, and a CS acquisition with an acceleration factor of R=12.4 − 13.8, using Cartesian pseudo-radial golden angle sampling pattern 39 and data driven-respiratory motion detection, as in 11 . Only data in expiration were kept for reconstruction as shown in Fig. 1.
To evaluate reconstruction accuracy on pathological anatomy, 4D flow data was acquired in a single patient with dilation of the ascending aorta, combined aortic stenosis and regurgitation due to a bicuspid aortic valve on a 3T Philips Ingenia system (Philips Healthcare, Best, the Netherlands) using a navigator-gated 2-fold accelerated parallel imaging 5 scan.
A receiver coil with 28 channels was used for acquisition which were reduced to 5 channels using coil compression 40 . Coil sensitivity maps were estimated with ESPIRiT 41 . Concomitant field correction was applied to the signal phase according to Bernstein et al. 42 and eddy currents were corrected for with a third-order polynomial model fitted to stationary tissue 43,44 .

Evaluation
We compared the proposed FlowVN to the state-of-the-art compressed sensing LLR-regularized (8) reconstruction 14 and the variational network by Hammernik et al. 15 which we refer to as HamVN. The LLR implementation from the Berkeley advanced reconstruction toolbox (BART) 45 was used with patch size p=8 and a maximum number of optimization iterations of 80. The optimal value of regularization parameter λ LLR =2.06 was chosen via grid search to minimize the reconstructed flow field residual v − v * 2 averaged over the manually segmented aorta on the retrospectively 12× undersampled acquisition. Since the original VN 15 was proposed for magnitude reconstruction of 2D and 3D data, we introduced the following modifications for the presented 4D flow evaluation: (i) 3D filters were grouped into 4 banks as in FlowVN (see Supplemental Algorithm 1), (ii) 1 -norm of reconstruction was optimized (i.e. equation (3) with τ →+∞ ). We refer to this architecture as "HamVN", the number of network layers, filters and control knots were the same as in FlowVN.

Retrospective Study
For simulated retrospective undersampling experiments, we used 2× PI data and simulate a pseudoradial Golden angle sampling pattern 39 with acceleration factors of 6 to 22.
For each undersampling factor we evaluated the nRMSE of the image magnitude, the relative error (RelErr) of velocity magnitudes inside the aorta and the angular error (AngErr) of the estimated velocity vectors: nRMSE(a, a * ) = N i (a i − a i ) 2 N max j (a j ) 2 , RelErr(a, a * ) = a − a * 2 a * 2 , Additionally, we report the structural similarity index (SSIM) 46 with σ SSIM =1.5 on the reconstructed magnitude images.

Prospective Study
Using manual aorta segmentations we compute flow over cross sections of the aorta by integrating velocity components projected onto the cross section normal. The peak flow is then defined as the maximal flow over cardiac phases for a given cross section. Moreover, we calculate peak throughplane velocity defined as maximum velocity projection across cross sections of the aorta over cardiac phases.
To quantify agreement with the reference 2× PI reconstruction, we performed Bland-Altman analysis 47 of peak flow and peak through-plane velocities.

Data and code availability
The code for the network training and inference used in this study and network weights are avaliable on CodeOcean together with volunteer data: https://codeocean.com/capsule/0115983/tree 48 Figure 6: Retrospective ablation study. a, Comparison of VN architectures and modifications proposed in this paper, using retrospectively undersampled data from 7 healthy volunteers for different acceleration factors R. HamVN-3D performs 3D filtering of yzt dimensional data (assuming that data is fully sampled in the x dimension). Modifications are accumulated as presented in the legend and, eventually, yield the FlowVN configuration. Undersampling-dependent modulation improves reconstruction quality only for higher acceleration factors. Note that the layer-wise exponential weighting of the reconstruction loss does not increase model complexity and, therefore, does not improve reconstruction accuracy for retrospective acceleration simulations. b, Velocity magnitude error maps for R=16 shown at systolic peak flow.