To develop a physics‐guided deep learning (PG‐DL) reconstruction strategy based on a signal intensity informed multi‐coil (SIIM) encoding operator for highly‐accelerated simultaneous multislice (SMS) myocardial perfusion cardiac MRI (CMR).
First‐pass perfusion CMR acquires highly‐accelerated images with dynamically varying signal intensity/SNR following the administration of a gadolinium‐based contrast agent. Thus, using PG‐DL reconstruction with a conventional multi‐coil encoding operator leads to analogous signal intensity variations across different time‐frames at the network output, creating difficulties in generalization for varying SNR levels. We propose to use a SIIM encoding operator to capture the signal intensity/SNR variations across time‐frames in a reformulated encoding operator. This leads to a more uniform/flat contrast at the output of the PG‐DL network, facilitating generalizability across time‐frames. PG‐DL reconstruction with the proposed SIIM encoding operator is compared to PG‐DL with conventional encoding operator, split slice‐GRAPPA, locally low‐rank (LLR) regularized reconstruction, low‐rank plus sparse (L+S) reconstruction, and regularized ROCK‐SPIRiT.
Results on highly accelerated free‐breathing first pass myocardial perfusion CMR at three‐fold SMS and four‐fold in‐plane acceleration show that the proposed method improves upon the reconstruction methods use for comparison. Substantial noise reduction is achieved compared to split slice‐GRAPPA, and aliasing artifacts reduction compared to LLR regularized reconstruction, L+S reconstruction and PG‐DL with conventional encoding. Furthermore, a qualitative reader study indicated that proposed method outperformed all methods.
PG‐DL reconstruction with the proposed SIIM encoding operator improves generalization across different time‐frames /SNRs in highly accelerated perfusion CMR.
Myocardial perfusion cardiac MRI (CMR) is used for functional assessment of stenoses in diagnosing coronary artery disease (1)(2)(3)(4)(5)(6)(7). Clinically, myocardial perfusion CMR is acquired using snap-shot imaging during the first pass of an exogenous contrast agent, which results in limited resolution and coverage (8)(9)(10). Low spatial resolution has been associated with dark rim artifacts that can compromise assessment of perfusion abnormalities (11). Additionally, coverage is typically limited to 3-to-4 noncontiguous slices (12), which may result in missed regions in microvascular disease. Furthermore, limited temporal resolution is associated with low contrastto-noise ratios and may produce cardiac motion artifacts (13). Therefore, trade-offs between spatio-temporal resolution and coverage still remain a major challenge in myocardial perfusion CMR, necessitating accelerated imaging techniques.
Parallel imaging has long been utilized in perfusion CMR but is limited to 2-to-3-fold acceleration (12). Spatio-temporal reconstruction (k-t) methods (14)(15)(16) have been proposed, but their acceleration rates remained limited (17). Subsequently, compressed sensing, low-rank methods, and their combinations have been adopted to perfusion CMR reconstruction to enable higher acceleration rates . These have enabled 3D whole heart myocardial perfusion (41)(42)(43)(44)(45)(46)(47)(48)(49), although a recent study has shown that 2D high resolution scans with smaller temporal footprint are more sensitive for detecting ischemia (50). Recently, simultaneous multislice (SMS) imaging has gained interest in CMR for improved coverage with minimal loss in image quality and signalto-noise ratio (SNR) (21,23,(51)(52)(53)(54). Yet, ultra-high acceleration rates are still limited when SMS imaging is combined with in-plane acceleration due to noise amplification (55).
Physics-guided deep learning (PG-DL) techniques have recently gained substantial interest in accelerated MRI, showing improved reconstruction quality at high acceleration rates compared to parallel imaging or compressed sensing (56)(57)(58)(59)(60)(61)(62)(63). These PG-DL techniques use a forward encoding operator incorporating MRI physics, while the proximal operation associated with regularization is solved implicitly by neural networks (61). However, PG-DL networks have several challenges that hamper their applicability in perfusion CMR. A 2D implementation processing slices/timeframes individually is a natural choice from an implementation perspective, and for avoiding temporal blurring. However, signal intensity changes across time-frames hinder the utility of such PG-DL networks, which have exhibited generalizability issues with such variations (64). An alternative way to train PG-DL reconstruction for perfusion CMR would be using a spatiotemporal network, yet this has its own challenges including memory limitations (65) and difficulty of procuring high-quality training databases due to differences in contrast uptakes/breathing patterns among subjects. Thus, application of PG-DL reconstruction to perfusion CMR has been difficult, and existing DL methods for perfusion CMR reconstruction have been limited to datadriven image enhancement networks (66)(67)(68), which are trained in a supervised manner using conventional compressed sensing reconstruction outputs as reference images. While this line of work improves reconstruction speed, the reconstruction quality is inherently limited by the conventional reconstruction used as reference for supervised training, which in turn hinders the true potential of DL reconstruction for perfusion CMR.
In this study, we propose to use a signal intensity informed multi-coil (SIIM) encoding operator in PG-DL networks to improve highly-accelerated perfusion CMR reconstruction. The proposed SIIM encoding operator is inherently aware of contrast/SNR changes across time-frames, leading to a uniform/flat signal level at the output of the network, which in turn assists the generalizability of PG-DL methods. Proposed SIIM encoding operator was compared with PG-DL using conventional operator, and conventional reconstruction methods, including split slice-GRAPPA (69), locally low-rank (LLR) regularization (34,70), regularized ROCK-SPIRIT (71) and a lowrank plus sparse (L+S) reconstruction (35) for free-breathing first-pass perfusion with 3-fold SMS and 4-fold in-plane acceleration. Results show that PG-DL reconstruction with the proposed SIIM encoding operator improves upon the other methods by reducing noise and residual artifacts.
The inverse problem for MRI reconstruction is formulated as an optimization problem
where 𝐲𝐲 Ω is the acquired multi-channel k-space, Ω is the in-plane undersampling pattern, 𝐄𝐄 Ω is the multi-coil encoding operator, 𝐱𝐱 is the image of interest, and 𝐧𝐧 is measurement noise. At high acceleration rates, this system is typically ill-conditioned. The first quadratic term enforces the data fidelity (DF) with acquired k-space points, and the second term ℛ(⋅) is a regularizer. This objective function may be solved using a multitude of techniques (72), which decouple the DF and regularizer terms into a series of sub-problems, including variable splitting with quadratic penalty (61), described in detail in Supporting Figure S1..
The encoding operator 𝐄𝐄 Ω in [1] is given as:
where 𝐅𝐅 Ω is a sub-sampled Fourier operator sampling the k-space locations specified by Ω, and 𝐒𝐒 c is a diagonal matrix representing the 𝑐𝑐 𝑡𝑡ℎ coil sensitivity map. In practice, 𝐒𝐒 c are estimated via ESPIRiT (73), and inherently encode B 1 -, which remain fixed across time-frames. Therefore, the solution of [1] presents varying signal intensities across time-frames, which mirror SNR variations in acquired k-space across time-frames.
We propose to encode dynamically-varying signal intensity in the encoding operator for PG-DL reconstruction. Let 𝐋𝐋 be a diagonal matrix whose entries are the pixel values of an image that contains the signal intensity information of a given time-frame. We define the SIIM encoding operator as:
where the inherent signal intensity variation across time-frames is encoded into encoding operator via 𝐋𝐋. Note that for perfusion CMR, we indeed have multiple 𝐋𝐋 𝐭𝐭 , t ∈ {1, ⋯ , T } where T is the number of time-frames, but for ease of notation, we use 𝐋𝐋 for a given time-frame of interest.
Consequently, the inverse problem for SIIM encoding operator is:
In the absence of a regularizer, it is easy to show (74)
where * is the Hermitian transpose. Thus, the underlying signal intensity information is restored by multiplication with the corresponding signal intensity informed images.
Signal intensity variations for a given time-frame can be captured with a low-resolution image, generated from central k-space, as the diagonal entries of 𝐋𝐋. In the context of parallel imaging, a similar concept was utilized, where low-resolution images from central k-space were used as coil maps, without normalizing them by their root-sum-squares image (74), and the signal intensity information was restored by multiplication with the root-sum-squares image, as in Eq. [5]. In this work, we instead use the formulation in Eq. [3], since this enables a more synergistic combination with ESPIRiT map estimation.
There are two major differences between SIIM and conventional encoding operators. First, there are numerical differences in solving the objective functions in [1] and [4], which was also noted for the unregularized case in parallel imaging (74). Thus, the SIIM formulation may overcome numerical instabilities at high acceleration rates. Second, in the regularized setup, the SIIM encoding operator has the additional benefit that the solutions of [4]
Free-breathing first-pass myocardial perfusion CMR was acquired on a 3T Siemens Magnetom Prisma (Siemens Healthineers, Erlangen, Germany) in 8 subjects (6 men, 2 women, age:39±18 years). This study was approved by our institutional review board, and written informed consent was obtained before each examination. A saturation-prepared GRE sequence was used, with relevant imaging parameters: FOV=360×320mm 2 ; spatial resolution=1.7×1.7mm 2 ; slice thickness=8mm; temporal resolution=116ms; SMS factor=3 (1/3 FOV shifts with CAIPIRINHA ( 75)); in-plane acceleration=4 (uniform undersampling, no ACS) and partial Fourier=6/8 (overall 16-fold acceleration) (55). Non-prepared GRE calibration scans were acquired at a lower spatial resolution=1.7×5.6mm 2 individually for all 9 slices. Details of the imaging sequence are given in Supporting Table S1.
Coil maps (𝐒𝐒 C ) were generated via ESPIRiT using central 24×24 regions of the calibrations scans of the corresponding slices (73) . Low-resolution images (𝐋𝐋) for each time-frame and slice were generated from the central 24×24 k-space region reconstructed using split slice-GRAPPA (69).
Note that this intermediate reconstruction step was necessary due to the lack of individual k-spaces for the slices of individual time-frames resulting from SMS encoding, and would not be necessary for single-slice/volume imaging. Subsequently, a Blackman filter was applied for ringing (76), followed by taking the magnitude of the SENSE-1 combination of individual coil images (77).
Finally, SIIM encoding operator 𝐇𝐇 Ω was generated by multiplying 𝐄𝐄 Ω by 𝐋𝐋, whose diagonal entries were the intensity values of the aforementioned magnitude SENSE-1 image, as in Eq. [3]. Further implementation details for SMS encoding are provided in Supporting Figure S2.
Due to lack of fully-sampled reference data in this highly-accelerated SMS perfusion CMR acquisition, the recently proposed self-supervised learning via data undersampling (SSDU) was used for training (61,78,79). Details of the multi-mask version of SSDU (80) S3.
Additionally, a numerical perfusion phantom (81) was used to evaluate the performance of different reconstruction methods, using in-vivo trained models. The details and results of these numerical experiments are presented in Supporting Tables S4-5 and Figures S5-7.
Qualitative image quality assessment was performed by an experienced cardiologist (15 years of experience). The reader was blinded to the reconstruction methods, orders of which were randomized. 4 Difference images between various reconstructions and linear baseline reconstruction split slice-GRAPPA are depicted in Supporting Figure S8. Proposed method shows noise-like differences with respect to split slice-GRAPPA, whereas residual artifacts are seen in all other regularized reconstructions. Videos of two subjects are included in Supporting Video S1-S2. Similarly, the proposed method shows the least amount of blurring and highest perceived SNR among all methods.
In this study, we proposed SIIM encoding operator for PG-DL reconstruction of image series with varying contrast across time-frames, and applied it to highly-accelerated myocardial perfusion CMR. The main advantage of using SIIM encoding operator is a uniform/flat signal level across different time-frames at the unrolled neural network output. This in turn facilitates generalizability of PG-DL reconstruction. The proposed approach improved upon multiple regularized reconstructions, showing better image quality, and reduced noise amplification and aliasing.
Conventional and SIIM encoding have two main differences. First, the solution of [1] using conventional encoding is adversely affected by ill-conditioning at high accelerations (83). As noted earlier, a similar concept to SIIM encoding was proposed in (74) to improve numerical stability for parallel imaging. The proposed SIIM encoding in [3] aims for a similar improvement, while enabling a synergistic combination with ESPIRiT, thus not necessitating a different coil map generation process as in (74). Second, and more importantly, in the PG-DL setup, SIIM encoding provides a more uniform contrast at the neural network outputs across time-frames. The output signal intensity of PG-DL with conventional encoding operator fluctuates across time-frames, and the regularization in the unrolled network needs to work with dramatically different signal levels.
On the other hand, SIIM encoding operator maintains a uniform output in terms of signal level.
Thus, regularization operates on more uniform SNRs in image space, for the corresponding outputs 𝐱𝐱 � SIIM , which empirically generalizes better across time-frames. Even though this intermediate solution has more uniform signal intensity, the final reconstruction is generated by multiplying with the corresponding low-resolution image for that time-frame, restoring the original signal intensity, as indicated in [2]. Thus, the use of SIIM encoding operator should not affect quantification in myocardial perfusion, consistent with conclusions from the uptake curves.
Finally, on a first look, Eq. [5] may resemble preconditioners in other MRI reconstruction problems (84)(85)(86), which are used to reduce the number of iterations for data fidelity. However, such preconditioners do not change the output signal intensity, thus solution of the preconditioned system coincides with that of the objective in [1]. Hence, SIIM operator is distinct from this typical use of preconditioning, leading to a more uniform signal intensity across time-frames.
DL reconstruction has gained interest in perfusion CMR, but has been limited to data-driven image enhancement approaches that learn a mapping between aliased and artifact-free images (66-68).
PG-DL approaches, which have been shown to outperform image enhancement methods (62,87,88) have remained elusive for perfusion CMR. One of the main challenges for PG-DL techniques has been related to generalizability with SNR changes (64), limiting the use of such reconstructions across perfusion time-frames, which is the main issue tackled in this study.
Another challenge for DL reconstruction in perfusion CMR has been the lack of gold-standard reference data. Aforementioned data-driven DL methods (66-68) were trained using supervision with compressed sensing reconstructions, limiting the performance of DL reconstruction. On the other hand, PG-DL methods enable self-supervised training from undersampled k-space data only (61,78,80,89,90), without a reference image. Thus, the combination of SIIM encoding and selfsupervised learning for PG-DL, as in this study, has the potential to further improve the utility of DL reconstruction for perfusion CMR (61). Finally, PG-DL reconstruction can be trained with fewer datasets compared to data-driven DL methods, and the number of k-spaces used for training in this work was in line with earlier PG-DL works that used ~200-to-360 k-spaces (56)(57)(58)(59)(60)(61)64,91,92), and was gathered using only 4 subjects. We note that the performance gap between the DL methods may change with a substantially larger training database, but this could not be investigated with our current cohort size.
The use of SMS encoding in this study required several design choices related to calibration data.
First, since central k-space data was not available for individual slices for each time-frame, an initial split slice-GRAPPA reconstruction was used to generate L, which suppresses aliasing but shows noise amplification. However, since only a limited central k-space region, containing high-SNR low-frequency k-space points, was used to generate L, SNR reduction effects from split slice-GRAPPA were observed to be minimal in subsequent processing. We emphasize that this step was only needed because of SMS encoding, and is not necessary for conventional 2D/3D encoding, where central k-space can be fully-sampled. Furthermore, Blackman filtering was used to avoid ringing, and the reader study did not report any dark rim artifacts associated with the use of L.
Second, calibration data for SMS reconstruction was acquired separately in free-breathing, which may be in different respiratory/cardiac motion states than perfusion data. Previously, it was shown that there were no differences between using free-breathing and breath-held calibration in another SMS CMR application in healthy cohorts (93). Furthermore, ESPIRiT uses only a 24×24 central region, leading to smooth maps, where motion-related artifacts in coil estimation may be nonsevere for most cohorts. However, evaluation of these pre-acquired calibration scans warrants further investigation, especially in patient populations with pharmacologically induced stress.
Finally, uniform undersampling was used in combination with SMS, since it allows easier integration in clinical sequences, and enables comparisons with clinically-used split slice-GRAPPA reconstruction. We note that compressed sensing methods are often used with random undersampling, thus their performance with uniform undersampling may be deteriorated.
This study has several limitations. All acquisitions in this study were prospectively accelerated.
Therefore, there is no gold-standard reference for image quality assessment. Since it is difficult to acquire first-pass perfusion on subjects multiple times due to need for repeated contrast injection, a conventional low-resolution perfusion scan with limited coverage was not performed, excluding a more typical clinical baseline for comparison. Additionally, no stress imaging data was available, which is clinically imperative for perfusion diagnostics. A pixel-wise mapping of myocardial blood flow (MBF) estimation ( 94) may be performed for quantitative assessment ( 95), but such analyses typically require modifications to the imaging protocol, such as administering dual doses (96,97) or using dual sequences (98)(99)(100)(101). Thus, MBF estimation could not be reliably performed with our acquisition protocol. Further clinical studies are warranted to assess full potential of the proposed method, and its diagnostic value in patients with suspected coronary artery disease.
The proposed PG-DL reconstruction with SIIM encoding operator generalizes well across timeframes/SNRs, and substantially improves upon several existing reconstruction methods for highlyaccelerated perfusion CMR. The product of the network output for the SIIM operator with low resolution images (𝐋𝐋) yields similar contrast to the network output for the conventional operator. S1: Free-breathing first-pass myocardial perfusion CMR imaging sequence details (4). Three sets of SMS-accelerated slices were acquired for a total of 9 slices, covering the whole heart where OVS modules were interleaved between every 9 imaging pulses to maintain suppression throughout the imaging (4). A non-prepared GRE was used to acquire calibration scans with FOV = 360 × 360 mm 2 and at a lower spatial resolution = 1.7 × 5.6 mm 2 individually for all 9 slices in free-breathing.
Conjugate Gradient (20 iterations)
ResNet (Supporting Information Figure 4)
ESPIRiT Calibration Kernel 6×6
ESPIRiT Threshold 0.02 / �|Ω|/(n RO ⋅ n PE )�, where 𝑛𝑛 𝑅𝑅𝑅𝑅 and 𝑛𝑛 𝑃𝑃𝑃𝑃 being the image sizes, and 𝐶𝐶 the number of channels. Note that this preconditioner amounts to a constant multiple of the identity matrix in the conventional encoding operator setting, reducing to the conventional CG algorithm.
Coil Distance 350 mm (distance of coil centers from origin)
In Supporting Information Table S4: Implementation details of the numerical phantom. Other imaging parameters were set to match with the in-vivo experiments as in Supporting Information Table S1. For the reconstruction experiments, the in vivo trained models were used for both PG-DL methods with conventional and SIIM encoding, which also test their generalizability to a phantom with different image features than in vivo images. Thus, no additional trainings were performed. Rest of the reconstruction techniques were performed with the details given in Supporting Information Table S3. Peak signal-tonoise ratio (PSNR) and structural similarity index measure (SSIM) were calculated between single-band fully sampled reference images and all reconstruction techniques. Error images were calculated as the difference between reference images and the reconstructed images.