Residual RAKI: A hybrid linear and non-linear approach for scan-specific k-space deep learning.

Parallel imaging is the most clinically used acceleration technique for magnetic resonance imaging (MRI) in part due to its easy inclusion into routine acquisitions. In k-space based parallel imaging reconstruction, sub-sampled k-space data are interpolated using linear convolutions. At high acceleration rates these methods have inherent noise amplification and reduced image quality. On the other hand, non-linear deep learning methods provide improved image quality at high acceleration, but the availability of training databases for different scans, as well as their interpretability hinder their adaptation. In this work, we present an extension of Robust Artificial-neural-networks for k-space Interpolation (RAKI), called residual-RAKI (rRAKI), which achieves scan-specific machine learning reconstruction using a hybrid linear and non-linear methodology. In rRAKI, non-linear CNNs are trained jointly with a linear convolution implemented via a skip connection. In effect, the linear part provides a baseline reconstruction, while the non-linear CNN that runs in parallel provides further reduction of artifacts and noise arising from the linear part. The explicit split between the linear and non-linear aspects of the reconstruction also help improve interpretability compared to purely non-linear methods. Experiments were conducted on the publicly available fastMRI datasets, as well as high-resolution anatomical imaging, comparing GRAPPA and its variants, compressed sensing, RAKI, Scan Specific Artifact Reduction in K-space (SPARK) and the proposed rRAKI. Additionally, highly-accelerated simultaneous multi-slice (SMS) functional MRI reconstructions were also performed, where the proposed rRAKI was compred to Read-out SENSE-GRAPPA and RAKI. Our results show that the proposed rRAKI method substantially improves the image quality compared to conventional parallel imaging, and offers sharper images compared to SPARK and ℓ1-SPIRiT. Furthermore, rRAKI shows improved preservation of time-varying dynamics compared to both parallel imaging and RAKI in highly-accelerated SMS fMRI.

Introduction

Magnetic resonance imaging (MRI) remains one of the most important modalities for neuroimaging (Bandettini, 2012; Leung et al., 2015), but it still faces challenges due to lengthy scan times. Parallel imaging (PI) is frequently used in clinical MRI to reduce scan times. These methods utilize differences in the profiles of multiple receiver coils to reconstruct undersampled data (Griswold et al., 2002; Pruessmann et al., 1999; Sodickson and Manning, 1997). However PI methods inherently exhibit trade-offs between acceleration rate and noise amplification (Aja-Fernández et al., 2014).

Recently, machine learning techniques have received substantial interest for MRI reconstruction. Several of these methods build on the regularized inverse problem from compressed sensing (Lustig et al., 2007), and learn a machine learning-based regularizer (Aggarwal et al., 2018; Chen et al., 2018; Dar et al., 2020a, 2020b; Eo et al., 2018; Hammernik et al., 2018; Han et al., 2018; Hosseini et al., 2020a; Mardani et al., 2018; Qin et al., 2018; Quan et al., 2018; Schlemper et al., 2018; Wang et al., 2016; Yaman et al., 2020a; Yang et al., 2018, 2016). Other data-driven approaches find a direct mapping from undersampled data into a de-aliased image (Han et al., 2020; Lee et al., 2018; Zhu et al., 2018). However, most of these methods rely on large training databases of fully-sampled data for training. Training on data independent from reconstruction target may lead to potential risks, such as bringing unwarranted information into reconstruction result or losing fine details, especially if the training databases do not include sufficiently many examples with pathologies of interest (Knoll et al., 2019, 2020a; Muckley et al., 2020).

An alternate line of work considers scan-specific training. The first method of this nature, Robust Artificial-neural-networks for k-space Interpolation (RAKI) uses convolutional neural networks (CNNs) trained from ACS data for scan-specific k-space interpolation (Akçakaya et al., 2019). RAKI showed improvements in image quality and noise reduction compared to conventional PI methods using linear k-space interpolation. The scan-specific design allows RAKI to be employed when massive fully sampled training data are not available, such as coronary MRI (Hosseini et al., 2020b) and highly accelerated simultaneous multi-slice/multi-band (SMS/MB) MRI (Nencka et al., 2021; Zhang et al., 2019a, 2018a). It has also been applied to algorithms that rely on completion of locally low-rank k-space neighborhoods (Kim et al., 2019), and in approaches that utilize the whole sub-sampled k-space for reconstruction (Zhang et al., 2019d). However, while the reconstruction improvement is visually easy to identify in these methods, the source of improvement from RAKI-type approaches is difficult to interpret. With the increasing importance of interpretability in machine learning (Reyes et al., 2020; Vellido, 2020), it is important to develop explainable methods that will aid in clinical translation for the RAKI framework.

In this paper, we propose residual RAKI (rRAKI) as a machine learning approach with improved interpretability for k-space interpolation, combining the advantages of both linear PI and nonlinear machine learning reconstructions. It uses a residual CNN (ResNet) architecture (He et al., 2016) with linear convolutions on the skip connection. The linear skip connection implements a linear reconstruction similar to GRAPPA, creating a baseline reconstruction. The multi-layered CNN compensates for the imperfections arising from this linear component, such as noise amplification and residual artifacts. Both components are trained on scan-specific ACS data concurrently. The performance of proposed rRAKI, is studied in T2 -weighted and FLAIR brain imaging from fastMRI database (Knoll et al., 2020b), high-resolution 3T and 7T anatomical imaging, as well as highly accelerated SMS/MB fMRI at 3T. Our results show rRAKI has noticeable advantages in terms of both noise resilience and artifact removal over linear methods, such as GRAPPA, as well as better interpretability and sharper images compared to RAKI.

Materials and methods

Linear k-space interpolation

GRAPPA is one of the most commonly used linear PI reconstruction approaches for multi-channel MRI reconstruction. A set of linear shift-invariant convolution kernels are calibrated from Nyquist-rate sampled auto calibration signal (ACS) data, to interpolate skipped data from adjacent sampled lines. Let s (k k i) denote an acquired data point in k-space location (k k of the ith coil, R be the acceleration rate, n be the number of coils, and for notational convenience, let denote a neighborhood around k-space location (k k across all coils that includes the sampled points for integers D and D specifying the pre-determined kernel size. Let be a column vector whose entries are the corresponding k-space values of s(k k i) for every point in . Under this notation, GRAPPA reconstruction can be expressed as: where g are the linear convolution weights for estimating the mth skipped line of coil i, for 1 ≤ m ≤ R-1. These convolution kernels, g are estimated by solving a least squares problem prior to reconstruction from the ACS region in k-space, where the ACS data are utilized as both regressor and regresand in Eq. (2)

Nonlinear RAKI reconstruction

RAKI achieves scan-specific nonlinear k-space interpolation by replacing the linear convolution kernels with several compact CNNs, each consisting of convolutions and nonlinear activations. CNNs in RAKI are trained from ACS data in the similar fashion of GRAPPA. For processing, the complex k-space is mapped to the real field, expanding the data to 2n input channels in total. In RAKI, each channel has its assigned CNN that generates R-1 k-space lines between sampled lines, while other possibilities of input-output designs can still be considered for reduced computational load (Zhang et al., 2019b). Let be the same neighborhood as the one in Eq. (1), only now defined over 2n real channels. Let be the R-1 missing lines between two sampled positions, adjacent to (k k) in channel j. Similarly, we define to be a column vector with its entries are elements s (α, β, γ), where . RAKI estimates the missing lines in k-space from acquired lines using CNNs by: where f (•) represents a CNN estimating the unacquired points in channel j based on the acquired data from all channels. f (•) was designed as a three-layered CNN in (Akçakaya et al., 2019). In practice, the CNN architecture can be designed in various forms in an application-specific manner (Zhang et al., 2019b). Point-wise nonlinearity, ReLU(x) = max(x, 0), is applied to all convolutional layers before the output layer. The last layer only conducts convolution to generate the final estimate.

Proposed rRAKI reconstruction

Fig. 1 depicts the residual network architecture used in rRAKI. Both the linear convolution component and the CNN take the sampled lines as input. The outputs from both components are combined into the final output, giving an estimation of all the skipped lines for a given channel, similar to RAKI. Let G denote the linear convolution component and F denote the nonlinear CNN component of the rRAKI network for channel j. Using this notation, the rRAKI reconstruction for a given channel j is summarized as: G and F are trained using ACS data, similar to (Akçakaya et al., 2019). rRAKI is designed so that G captures the linear reconstruction baseline, while F reduces the residual artifacts and noise amplification non-linearly. Thus, the training, which is performed over the ACS region, aims to minimize the error for the estimation that combines the G and F components, while also minimizing the error associated with the G component. Let y be the vector containing the elements of corresponding to the target points in the ACS region of channel j, where and specify the dimension of ACS along x and y axis, and let y be the vector containing the elements of corresponding to the source points across all channels in the ACS region. Training of the rRAKI network is performed using the following loss function: where || • ||2 is the l2 norm, and are the trainable weights in the G and F networks respectively, and λ is a weighting factor. Further implementation details are provided in Section 2.7. The reconstruction is then performed using the learned parameters and

Fig. 1.

Schematic of the rRAKI network. The network takes undersampled k-space data from all coils as input, and interpolates the k-space of a certain coil. It consists of two branches: the linear part G, which provides a linear reconstruction baseline, and the nonlinear part F, which achieves nonlinear noise and artifacts removal from G. G can be either a single convolution or a linear convolutional network. F is a nonlinear CNN that consists of linear convolutions and nonlinear activations. rRAKI network outputs G + F as the final reconstruction of the given channel. The outputs from G and F have R -1 channels corresponding to R -1 skipped lines between two sampled phase-encoding positions.

fastMRI brain datasets

Fully-sampled axial T2 (T2 -weighted) and FLAIR datasets from the fastMRI database (Knoll et al., 2020b) were employed. 300 slices of fully sampled data were retrospectively undersampled along phase-encode direction at rate 4. Reconstructions using GRAPPA, Tikhonov-regularized GRAPPA, non-linear GRAPPA (Chang et al., 2012), ℓ1 wavelet regularized SPIRiT (Lustig and Pauly, 2010), SPARK (Arefeen et al., 2022), RAKI and rRAKI were performed. For displaying and assessment purposes, multi-coil images were combined using root-of-sum-of-squares. SSIM, NRMSE, and blur metrics (Crété-Roffet et al., 2008) were computed with respect to the reference. Normality of SSIM, NRMSE and blur metrics were assessed using the Jarque-Bera test, prior to testing statistical differences in SSIM, NRMSE and blur metrics using paired t-tests. P-values < 0.05 were considered significant.

Anatomical imaging

Anatomical brain imaging was performed on a 3T Siemens Magnetom Prirma system and a 7T Siemens Magnex Scientific (Siemens Healthcare, Erlangen, Germany) system using 32-channel receiver head coil-arrays. The imaging protocols were approved by the local institutional review board, and written informed consent was obtained from all participants before each examination for this HIPAA-compliant study. For 3T imaging, a T1 -weighted 3D-MPRAGE sequence was acquired in a healthy subject with the following parameters: Field-of-view (FOV) = 224 × 224 × 179 mm3, resolution = 0.7 × 0.7 × 0.7 mm3, matrix size = 320 × 320, TR/TE = 2400 ms/2.2 ms, flip angle = 8°, bandwidth = 210 Hz/pixel, inversion time = 1000 ms, ACS lines = 40, with iPAT = 2 and 5. Furthermore, the R = 2 acquisition was also retrospectively undersampled to R = 4 and 6 (Akçakaya et al., 2019). For 7T imaging, 3D-MPRAGE was acquired in a healthy volunteer with the following parameters: FOV = 230 × 230 × 154 mm3, resolution = 0.6 × 0.6 × 0.6 mm3, TR/TE = 3100 ms/3.5 ms, flip angle = 6°, bandwidth = 140 Hz/pixel, inversion time = 1500 ms, ACS lines = 40, with R = 3, 4, 5, 6. Additionally, two averages were acquired for R = 5 and 6 data to mitigate the SNR loss from undersampling (Akçakaya et al., 2019). The k-space data were inverse Fourier transformed along the slice direction for all datasets, and a central slice was processed. Reconstructions were performed using GRAPPA, RAKI and rRAKI. All methods were calibrated or trained on the slice-specific ACS region that was built-in to the acquisition. For display and evaluation, multi-coil images were combined using root-of-sums-of-squares. Reconstruction quality was assessed qualitatively, since a fully-sampled reference was not available in these acquisitions.

Simultaneous multi-slice fMRI

SMS/MB fMRI data was acquired on a 3T Siemens Magnetom Prisma (Siemens Healthcare, Erlangen, Germany) scanner with a 32-channel receiver head coil-array. The Human Connectome Project protocol (Van Essen et al., 2012) was used with SMS/MB factor = 8 and blipped-CAIPI encoding (Setsompop et al., 2012) with a FOV/3 shift between adjacent multiband slices, resolution = 2 × 2 × 2 mm3 and TE/TR = 37/1000 ms, FOV = 208 × 180 × 144 mm3, flip angle = 52°, matrix size = 104 × 90, bandwidth = 2290 Hz/pixel. Calibration data containing the individual slices was acquired integrated with and prior to the fMRI image series at the same resolution. In order to test the potential of rRAKI at high SMS/MB acceleration rates, a SMS/MB = 16 acceleration was retrospectively simulated using this acquisition. Specifically, each fMRI acquisition contains nine slice groups of SMS/MB = 8 slices for a total of 72 slices, resulting in five retrospective SMS/MB = 16 slice groups for each subject. The methodology for simulating SMS/MB = 16 data from SMS/MB = 8 acquisition are further detailed in supplementary materials.

Reconstructions were performed using the proposed rRAKI, as well as linear parallel imaging via RO-SENSE-GRAPPA (RSG) (Moeller et al., 2010) and RAKI as comparison. Readout concatenation (Moeller et al., 2010) was used for all the reconstructions. In this method, for kernel calibration, individual unaccelerated images are concatenated in image domain along the readout direction, and then transformed into k-space via Fourier transform, where SMS/MB encoding can be viewed as acceleration in this concatenated readout direction. The concatenated unaccelerated images are then used for estimating the interpolation rule, e.g. linear convolutions in GRAPPA, in this extended space (Demirel et al., 2021; Moeller et al., 2010). In order to reduce overfitting issues for the GRAPPA-type convolution from a single calibration frame, additional calibration data was generated from the baseline SMS/MB = 8 reconstructions from the first 8 time frames. In supplementary materials we provide a detailed description of the calibration data generation for SMS/MB = 16 experiments in this study. For all tested methods, 8 calibration frames were used for generating the interpolation kernels or CNNs. Additionally, this amount of calibration data enables G in rRAKI to be implemented as a multi-layered linear convolutional network (Bell-Kligler et al., 2019), which keeps the linearity while enabling multiple optimal solutions and avoiding issues with local minima. Following the calibration stage, each individual time frame of the fMRI series were reconstructed using the same convolutional kernel and/or CNN. For display and evaluation, multi-coil images were combined using root-of-sums-of-squares. Due to the lack of fully-sampled data, reconstruction quality was assessed visually for individual frames. Further quantitative evaluation was performed using temporal SNR (TSNR), which measures the variation of image signal along time in a point-wise fashion. TSNR maps were calculated as the mean of the image series over time divided by its standard deviation (Tabelow et al., 2009).

Implementation details

GRAPPA and RSG were implemented using MATLAB 2016a (Math-Works, USA). RAKI and rRAKI was implemented using Tensorflow 1.7.0 and Python 3.6.2, supported by CUDA 8.0 and CuDNN 7.0.5, on Linux kernel 3.10.0. The Python environment was created under Anaconda 5.1.0. Adam (Kingma and Ba, 2015) was employed as the optimizer for network training. All programs were run on a server with two Intel E5–2643 CPUs (6 cores each, 3.7 GHz), 256 GB memory and two NVIDIA Tesla V100 GPU (32 GB memory each) with single precision.

For fastMRI dataset and anatomical imaging, the parameters of all non-machine learning reconstruction methods were empirically tuned for best performance. In addition to visual assessment, numerical metrics including SSIM and NRMSE were also considered in parameter tuning, when a reference image was available. GRAPPA and Tikhonov regularized GRAPPA was implemented using a 5 × 4 kernel (Griswold et al., 2002). An regularization factor of 0.05 was employed in Tikhonov regularized GRAPPA. ℓ1 -SPIRiT used a 5 × 5 kernel, and ℓ1 Daubechies-4 wavelet regularization. Regularized SPIRiT reconstruction was solved iteratively using ADMM. For the CNNs in RAKI and nonlinear part of rRAKI, we express the shape of convolution kernel used in layer 1 as , representing a convolution, where and denote the sizes of convolution window along k and k directions respectively, n and n denote the input and output channel number, respectively. Under this notation, n0 = 2n since the CNN takes k-space data of 2n channels. The output layer always outputs R – 1 channels for k-space estimation. In anatomical imaging, we employed the following hyper-parameters: w1 = [5, 2, 32], w2 = [1, 1, 8], w3 = [3, 2, R-1] for RAKI (Zhang et al., 2018a). rRAKI was implemented with a 5 × 2 linear convolution kernel, and its nonlinear part shared the same parameters as RAKI for all tested scenarios for comparison purposes. λ has been set to 1 for an equally weighting between linear and non-linear components. Parameters of Adam optimizer were set as α = 0.0003, β1 = 0.9, β2 = 0.999, ε = 10−8. SPARK was implemented as is suggested in (Arefeen et al., 2022).

For SMS/MB imaging, larger kernels were employed for RSG, RAKI and the proposed rRAKI due to the readout-concatenated nature of k-space. Specifically, RSG applied an 11 × 10 kernel. RAKI employed 5-layered CNNs with hyper-parameters: w1 = [11, 10, 32], w2 = [1, 1, 64], w3 = [3, 2, 32], w4 = [1, 1, 64], w5 = [5, 4, R-1]. rRAKI employed a multi-layered linear convolution network that has similar structure as F but without nonlinear activations. Both G and F shared the same hyper-parameters as those were used for RAKI, with λ = 1. Parameters of Adam optimizer were set as α = 0.0003, β1 = 0.9, β2 = 0.999, ε = 10−8.

Results

fastMRI dataset

Fig. 2 and Fig. 3 depict representative reconstruction results of T2 -weighted and FLAIR data from fastMRI database, respectively. ℓ1 -SPIRiT, RAKI, SPARK and rRAKI show visible advantages over GRAPPA, Tikhonov regularized GRAPPA and non-linear GRAPPA in terms of noise resilience. RAKI and rRAKI exhibit similar image quality, while providing visibly sharper images with fine details compared to ℓ1 -SPIRiT and SPARK. Note that ℓ1 -SPIRiT and SPARK exhibit more noise suppression albeit at the cost of smoother looking images. Table 1 lists the mean and standard deviation of SSIM, NRMSE and blur metrics for both T2-weighted and FLAIR datasets. For both T2-weighted and FLAIR datasets, all quantitative metrics showed normality, allowing subsequent paired t-tests to be conducted. The results of the paired t-tests are listed in supporting materials Table S1–S6 and summarized here. For both T2-weighted and FLAIR imaging, RAKI, rRAKI and SPARK show no statistical difference in terms of SSIM and NRMSE, while they are statistical different in terms of blur metrics, where SPARK had higher (worse) blur metrics than RAKI and rRAKI, indicating image blurring that matches the visual assessments. RAKI, rRAKI and SPARK outperform GRAPPA, Tikhonov-regularized GRAPPA, nonlinear GRAPPA and ℓ1-SPIRiT in terms of SSIM and NRMSE, with statistical significance. GRAPPA and Tikhonov-regularized GRAPPA show no statistical difference in all metrics in T2-weighted imaging, while for FLAIR imaging Tikhonov-regularized GRAPPA has statistically higher blur metrics value than GRAPPA. Similarly, nonlinear GRAPPA has statistically higher blur metrics than GRAPPA in T2-weighted imaging, while it has statistically better SSIM and NRMSE in FLAIR imaging. Nonlinear GRAPPA statistically outperforms Tikhonov-regularized GRAPPA in terms of NRMSE in T2-weighted imaging, While Tikhonov-regularized GRAPPA outperforms nonlinear GRAPPA in terms of SSIM and NRMSE in FLAIR imaging. ℓ1-SPIRiT outperforms GRAPPA, Tikhonov-regularized GRAPPA and nonlinear GRAPPA in terms of SSIM and NRMSE for both T2-weighted and FLAIR imaging, but ℓ1-SPIRiT also exhibits high blur metrics values indicating visual blurring.

Fig. 2.

Representative reconstruction results of T2-weighted images from the fastMRI database at R = 4, using GRAPPA, Tikhonov-regularized GRAPPA, non-linear GRAPPA, RAKI, SPARK and rRAKI. ℓ1-SPIRiT, RAKI, SPARK and rRAKI have reduced noise amplification compared to GRAPPA, Tikhonov-regularized GRAPPA and non-linear GRAPPA. ℓ1-SPIRiT and SPARK exhibit visual blurring and loss of fine details. RAKI and rRAKI exhibit visibly similar image quality, and provide sharper images and improved recovery of fine details compared to ℓ1-SPIRiT and SPARK.

Fig. 3.

Example results of FLAIR images from fastMRI database at R = 4, reconstructed using GRAPPA, Tikhonov-regularized GRAPPA, non-linear GRAPPA, RAKI, SPARK and rRAKI. ℓ1-SPIRiT, RAKI, SPARK and rRAKI outperform GRAPPA, Tikhonov-regularized GRAPPA and non-linear GRAPPA in terms of noise reduction. Blurring and loss of fine details are visible in ℓ1-SPIRiT and SPARK. RAKI and rRAKI demonstrate visibly similar image quality, while providing sharper images and improved recovery of fine details compared to ℓ1-SPIRiT and SPARK.

Table 1

Mean and standard deviations of SSIM, NRMSE and blur metrics of reconstruction methods tested on T2-weighted and FLAIR images from the fastMRI database at R = 4. Statistical analyses show that machine learning-based approaches RAKI, SPARK and rRAKI improve on ℓ1-SPIRiT GRAPPA and its variants for all metrics. Furthermore, RAKI and rRAKI statistically improve on SPARK in terms of the blur metric, matching visual observations. There is no statistical difference between RAKI, SPARK and rRAKI in terms of SSIM and NRMSE.

	T2-Weighted			FLAIR
	SSIM	NRMSE	Blur Metrics	SSIM	NRMSE	Blur Metrics
Reference	N/A	N/A	0.291 ± 0.0333	N/A	N/A	0.288 ± 0.0141
GRAPPA	0.868 ± 0.0611	0.110 ± 0.0241	0.244 ± 0.0283	0.901 ± 0.0299	0.0814 ± 0.0141	0.242 ± 0.0142
Tikhonov GRAPPA	0.873 ± 0.0455	0.111 ± 0.0211	0.248 ± 0.0273	0.902 ± 0.0267	0.0815 ± 0.0132	0.244 ± 0.0146
Non-linear GRAPPA	0.871 ± 0.0764	0.107 ± 0.0288	0.251 ± 0.0276	0.886 ± 0.0450	0.0875 ± 0.0201	0.243 ± 0.0170
ℓ1-SPIRIT	0.891 ± 0.0395	0.101 ± 0.0179	0.337 ± 0.0432	0.909 ± 0.0214	0.0766 ± 0.0135	0.319 ± 0.0176
RAKI	0.904 ± 0.0380	0.0904 ± 0.0161	0.278 ± 0.079	0.927 ± 0.0205	0.0674 ± 0.0121	0.271 ± 0.0134
SPARK	0.909 ± 0.0463	0.0893 ± 0.0168	0.305 ± 0.0334	0.926 ± 0.0223	0.0687 ± 0.0123	0.300 ± 0.0172
rRAKI	0.901 ± 0.0446	0.0905 ± 0.0168	0.269 ± 0.0257	0.928 ± 0.0221	0.0679 ± 0.0128	0.264 ± 0.0133

Computation times of GRAPPA, RAKI, SPARK and the proposed rRAKI are reportedin Supplementary Material Figure S3. The pure linear nature of GRAPPA allows a fast reconstruction speed, taking approximately a second to reconstruct the fastMRI brain image at rate 4. RAKI and rRAKI takes less than 20 s for the reconstruction. Owing to the use of a linear reconstruction baseline, rRAKI converges in fewer epochs than RAKI. Note as detailed in Section 2.7, GRAPPA was implemented in MATLAB, while the other methods were implemented using TensorFlow.

Fig. 4 depicts the results from the reconstruction of the 3T MPRAGE data using GRAPPA, RAKI, rRAKI, for different acceleration rates, as well as the linear part G and the nonlinear part F (scaled by 3 for improved visualization) of rRAKI. RAKI and rRAKI show minor advantages over GRAPPA in noise resilience below rate 4. The advantage becomes more noticeable for acceleration rates of 5 and 6, where RAKI and rRAKI have visually lower noise than GRAPPA. For rRAKI, the G part shows similar appearance to the GRAPPA results, in accordance with the interpretable nature of the reconstruction. The noise amplification in the linear G part is reduced using the nonlinear F part. The combination of these two components matches RAKI in terms of noise resilience.

Fig. 4.

Reconstruction results for a 3T MPRAGE acquisition using GRAPPA, RAKI, and rRAKI, for acceleration rates 4 to 6. For display purposes, the G component of rRAKI is filled with the acquired lines at the acquired positions, while the F component, which is the residual from G, has zeros in the acquired positions. At rate 4, RAKI and rRAKI perform similarly to GRAPPA in terms of noise resilience. This difference becomes more pronounced at higher rates 5 and 6. As expected, the linear portion (G) of rRAKI matches the image quality of GRAPPA. Amplified noise can be observed in G under all tested rates, which is subsequently reduced by F. Combined output from rRAKI matches the image quality of RAKI.

Similar observations apply to 7T MPRAGE data, the results of which are depicted in Fig. 5. For this dataset, all tested approaches successfully remove aliasing artifacts for the rates shown. The linear part of rRAKI and GRAPPA results present similar visual quality, followed by visible noise removal given by the nonlinear F part of rRAKI. rRAKI matches RAKI in terms of noise resilience, and both of them show visible advantage over GRAPPA, especially at rates 5 and 6, while rRAKI has improved interpretability.

Fig. 5.

Reconstruction results for a 7T MPRAGE acquisition using GRAPPA, RAKI and rRAKI for acceleration rates 4 to 6. G and F components of rRAKI are displayed as described in Fig. 4. At rate 4, there are no noticeable differences between the three methods. Visible noise amplifications are observed in GRAPPA at rates 5 and 6. At these rates, RAKI and rRAKI show improved noise resilience compared to GRAPPA. Although the linear part (G) of rRAKI suffers from similar noise amplification to GRAPPA, the non-linear part (F) successfully reduces the noise level, leading to an enhanced image quality that matches RAKI.

Additional reconstruction results of 3T and 7T anatomical imaging are shown in supplementary Fig. S4 and Fig. S5, including GRAPPA, Tikhonov-regularized GRAPPA, ℓ1-SPIRiT, RAKI, SPARK and rRAKI reconstructions at acceleration R = 6. Non-linear GRAPPA, ℓ1-SPIRiT, RAKI, SPARK and rRAKI show improved noise resilience compared to GRAPPA and Tikhonov-regularized GRAPPA. However, residual aliasing is visible in ℓ1-SPIRiT and SPARK. Non-linear GRAPPA, RAKI and rRAKI exhibit closer visual quality.

Representative results for high rate SMS/MB reconstruction are presented in Fig. 6. Six representative slices are displayed to demonstrate improved visualization. A display of all unaliased SMS/MB = 16 images are provided in Supplementary Fig. S6 and S7. Although no fully sampled reference exists in this setting, high-quality SMS/MB = 8 reconstruction results with RSG are provided for baseline comparison. At SMS/MB = 16, RSG displays visible artifacts and reconstruction noise, which is reduced in RAKI and rRAKI. rRAKI provides sharper results than RAKI, while reducing artifacts and noise amplification.

Fig. 6.

Reconstruction results from simulated SMS/MB = 16 data using RSG, RAKI and rRAKI are displayed for 4 representative slices out of 16. Due to the lack of fully sampled ground-truth data, the conventional high-quality SMS/MB = 8 reconstructions from RSG are employed as a baseline reference. G and F components of rRAKI are displayed as described in Fig. 4. At SMS/MB = 16, both RAKI and rRAKI show noticeable improvements compared to RSG, including lower noise and fewer reconstruction artifacts. However, RAKI suffers from visible blurring, which is ameliorated using rRAKI.

Representative TSNR maps are displayed in Fig. 7. For SMS/MB = 16, both RAKI and rRAKI demonstrate increased TSNR compared to RSG. However, RAKI provides a visibly homogeneous TSNR compared to the MB8 reference, which suggests a loss of sensitivity to temporal dynamics due to suboptimal generalizability. On the other hand, rRAKI shows slightly lower TSNR values compared to SMS/MB = 8 baseline, but the TSNR map from rRAKI preserves structural information.

Fig. 7.

TSNR maps corresponding to the slices in Fig. 5. RSG results at SMS/MB = 8 are used as baseline reference. At SMS/MB = 16, RAKI and rRAKI show higher TSNR values compared to RSG. However, RAKI leads to a homogeneous TSNR map, with values greater than the SMS/MB = 8 reference, indicating a risk of losing temporal dynamics. rRAKI outperforms RAKI by preserving more structural information in its TSNR map, while showing anatomically similar albeit lower TSNR values to the SMS/MB = 8 reference.

Discussion

In this study, we proposed rRAKI, which combines a conventional linear reconstruction with a neural network based nonlinear k-space reconstruction for improved image quality and interpretability. Undersampled k-space is interpolated using linear convolutions as a baseline, while the nonlinear CNNs further suppress the residual errors arising from this linear part. Such residual errors include g-factor noise amplification (Robson et al., 2008), and calibration errors due to noise in the regressor and regressand (Akçakaya et al., 2019; Chang et al., 2012). The proposed rRAKI has demonstrated noticeable improvement in image quality in high rate anatomical and SMS/ MB imaging.

An important contribution of rRAKI is the ability to interrogate the components in the reconstruction, where the linear component ensures that the unaliasing of signal is accurate, and the non-linear component captures inhomogeneous noise amplification and other artifacts. The explicit separation of these artifacts, as well as the inference of system conditioning are two features of the rRAKI algorithm that are useful in determining performance for high acceleration. The proposed rRAKI strategy exhibited improved image quality compared to conventional linear PI reconstruction in anatomical imaging at higher undersampling factors, as well as in highly accelerated SMS/MB imaging.

One of the advantages of RAKI-type methods have been its scan-specific nature, allowing them to be used in the absence of large training databases (Arefeen et al., 2022; Hosseini et al., 2020b; Kim et al., 2019; Nencka et al., 2021; Zhang et al., 2019a, 2018a). Recently, several other studies have aimed to develop alternative scan-specific deep learning methods. In (Yaman et al., 2021a), physics-guided DL reconstruction with algorithm unrolling was performed in a zero-shot manner building on self-supervised learning methods (Yaman et al., 2020b, 2021b). Another line of work that uses adversarial models relying on deep image prior (Ulyanov et al., 2018) has also been proposed (Korkmaz et al., 2021), showing interest in this growing area.

Several other works have also demonstrated noticeable gains by connecting linear and nonlinear methods (Arefeen et al., 2022; Dar et al., 2021). In the context of database training, (Dar et al., 2021) takes advantages of a hybrid linear and non-linear structure, leading to significant advantages over conventional physics-guided model in terms of artifacts and noise reduction. In the context of scan-specific approaches, SPARK trains a non-linear CNN minimizing the reconstruction error arising from a fixed baseline reconstruction (Arefeen et al., 2022). In one of our preliminary studies for this work (Zhang et al., 2019a), a fixed linear baseline was shown to lead to difficulty with convergence when fitting acquired k-space data to the reconstruction residual, due to varying signal intensity in k-space. Thus, in this work, we adapted a joint optimization of the linear and non-linear components for improved training. In (Arefeen et al., 2022), the training issue with a fixed baseline was resolved using a fine-tuned CNN and promising results. In this work, without imposing any condition on the CNN, we performed several comparisons between rRAKI and SPARK, and demonstrated the advantages of rRAKI in terms of both visual quality and numerical metrics, especially related to image sharpness and recovery of fine details.

The explicit separation of the linear and non-linear components in rRAKI, along with the joint training of these components, enables the training and testing to be performed on more similar signal energy levels compared to RAKI. The linear component fits most of the signal energy in the ACS data, while the non-linear CNN characterizes the residual. Using the shift-invariance of the linear component, this strategy extends to the outer k-space, while the non-linear nature of CNN is empirically important to adaptively fit such data. The non-linear component used in this study, consisting of ReLU after convolution without biases, can be expressed as a multiplication of matrices, where the ReLU activations lead to signal-dependent binary diagonal matrices (Ye et al., 2018). The signal dependency of this component empirically enables an adaptive trade-off in regularization without the need for an explicit weighting parameter as in earlier works for regularizing GRAPPA-type reconstructions (Weller et al., 2012), which is consistent with how CNNs are able to work across different noise levels in denoising problems (Gnanasambandam and Chan, 2020; Zhang et al., 2018b). Although the hybrid linear and non-linear design of rRAKI can lead to a low calibration error during training, over-fitting may occur when rRAKI is trained with limited ACS data.

For SMS/MB = 16, a single calibration frame only has 129,000 k-space points that can be used for training, whereas the network we employed has 730,120 parameters for each channel to be reconstructed. Thus, the additional calibration data improved both the linear and combined reconstruction of rRAKI, by enabling improved estimation of linear convolutions, as well as the use of a deep network.

For neuroimaging, mapping of resting-state networks into frequency-bands above 0.2 Hz with high spatial resolution is one of the needs for bridging the functional connectomics from fMRI and magnetic resonance encephalography (MREG), (Hennig et al., 2007). Efforts on interrogating these connections have been pursued with compressed sensing reconstructions for MB-EVI acquisition (Vakamudi et al., 2018) and fast sequences such as line-scanning for laminar fMRI (Yu et al., 2014), since conventional acquisition and undersampling strategies have been insufficient. rRAKI for fMRI demonstrate one potential approach to maintain conventional acquisition methods, while improving both the quality of the reconstruction and the temporal stability necessary for probing of such biological systems.

Our study has limitations. Regularization strength in Tikhonov regularized GRAPPA and the regularization parameters in ℓ1-SPIRiT were empirically tuned in this study. This included both visual assessment, and quantitative metrics such as SSIM and NRMSE when a reference image was available. For Tikhonov regularized GRAPPA, a high regularization strength suppressed noise amplification but potentially led to visible aliasing leftovers. Thus, in this case, we also included the reduction of aliasing artifacts in our visual evaluation instead of noise suppression alone. Unlike the fastMRI experiments, quantitative evaluation with respect to a reference were not provided for 3T and 7T anatomical imaging, as well as SMS fMRI, due to the lack of fully-sampled reference data. In this study, we focused on a scan-specific approach, thus the database learning methods were excluded from the comparison, For our implementations, hyperparameters of the CNNs were empirically adjusted, similar to previous studies on RAKI-type methods (Arefeen et al., 2022; Zhang et al., 2018a). Although we only focus on 1D Cartesian undersampling in this study, one may perform rRAKI on other uniform sampling patterns, such as 2D PE acceleration in 3D Cartesian sampling, by replacing Eq. (1) with a corresponding neighborhood selector.

Conclusion

In this study, we proposed a machine learning based k-space reconstruction approach, rRAKI, for enhanced image quality and improved interpretability. The efficacy of rRAKI in reducing noise and residual artifacts compared to conventional parallel imaging and RAKI was shown in both anatomical and functional imaging.