Achieving high-resolution 1H-MRSI of the human brain with compressed-sensing and low-rank reconstruction at 7 Tesla.

Low sensitivity MR techniques such as magnetic resonance spectroscopic imaging (MRSI) greatly benefit from the gain in signal-to-noise provided by ultra-high field MR. High-resolution and whole-slab brain MRSI remains however very challenging due to lengthy acquisition, low signal, lipid contamination and field inhomogeneity. In this study, we propose an acquisition-reconstruction scheme that combines 1H free-induction-decay (FID)-MRSI sequence, short TR acquisition, compressed sensing acceleration and low-rank modeling with total-generalized-variation constraint to achieve metabolite imaging in two and three dimensions at 7 Tesla. The resulting images and volumes reveal highly detailed distributions that are specific to each metabolite and follow the underlying brain anatomy. The MRSI method was validated in a high-resolution phantom containing fine metabolite structures, and in five healthy volunteers. This new application of compressed sensing acceleration paves the way for high-resolution MRSI in clinical setting with acquisition times of 5 min for 2D MRSI at 2.5 mm and of 20 min for 3D MRSI at 3.3 mm isotropic.

Introduction

Magnetic resonance spectroscopy (MRS) [1] has been one of the main motivations driving MR towards ultra-high field (⩾7T). With the recent advent of FDA approval and CE certification of 7T MR systems for clinical use, there is high interest in developing robust fast high-resolution MRS imaging (MRSI) methods to map the neurochemistry of human brain with greater anatomical details. MRS at 7T and higher field, benefits from increased sensitivity and spectral dispersion [2,3], which are the two intrinsic factors limiting the chemical information that can be extracted at the current clinical (⩽3T) fields. However, capitalizing the benefits of ultra-high field have proven more challenging for human MRS imaging (MRSI) compared to single voxel spectroscopy (SVS) [4]. Large B0 and B1 inhomogeneity over the human brain [5,6] results in nonuniform image quality across an extended field of view. Acquisition of the 4D (k,t) space by traditional clinical MRSI sequences [7] is slow due to long repetition time (TR) needed to sample adequately the time domain, and due to notably high specific-absorption rate (SAR) at ultra-high field which further limits the minimum TR. These problems prohibit the acquisition of spatial high-resolution MRSI to benefit from the higher signal-to-noise (SNR) at ultra-high field. Significant progress has been done in the recent years in accelerating high-resolution MRSI at ultra-high field [8] by employing sequences with low SAR, short TR and short echo time (TE). In particular, the pulse acquire (FID) excitation [9,10] has grown in popularity at ultra-high field due to the low flip angle excitation allowing for reduced TR without T1-weighting of the signal, it is not limited by the maximum voltage delivered by transmit coil in human, it has reduced SAR, and the very short TE reduces signal loss due to relaxation. The elimination of fat suppression pulses further reduces SAR and TR at ultra-high field but requires advanced reconstruction methods to decontaminate metabolic signal from skull lipid leakage [11-14]. The problem of MRSI acquisition time at ultra-high field was addressed with the use of spectrospatial encoding using echo-planar [15], concentric circles [16], rosettes [17] and spiral [18] trajectories, but also with parallel imaging methods [19,20], and more recently with a combination of non-Cartesian trajectories and parallel imaging [21].

Compressed sensing (CS) is a technique that allows for faithful recovery of an undersampled signal [22]. Its application to MRI enable acceleration of the acquisition by random k-space undersampling [23,24]. Although now being widely used in clinics [25,26] for standard MRI, CS has been scarcely explored for 1H-MRSI in humans, mostly at 3T [27-30] and in one study for single slice 2D MRSI at 9.4T [31]. However, 2D MRSI is typically performed with thick slices (8–10 mm) which lead to partial volume averaging, while isotropic high-resolution 3D coverage is preferable to match better the fine details of gray/white matter structure. Initial CS demonstrations of 1H-MRSI at 3T have been shown with volume localization [27-30] which excludes the scalp to avoid lipid artifacts, but also misses most of the lateral cortex. Compressed sensing has been successfully used to speed-up acquisition of X-nuclei MRSI with full field of view (FOV) coverage, such as 13C [32] and 31P [33], where spectra are more sparse, the large nuisance signals of fat and water are absent, and the effect of B0 inhomogeneity is reduced. On the other hand for whole brain CS 1H-MRSI at ultra-high field, proper separation between metabolites and dominant fat and water signals is of paramount importance [34,35]. Compared to parallel imaging, the random undersampling in CS has the advantage that it does not produce structured aliasing artefacts, is not affected by the g-factor penalty, does not require a calibration scan [23], but instead involves a more lengthy reconstruction which is acceptable when a real-time or fast answer is not needed such as often is the case for MRSI, although clinically acceptable reconstruction times can be obtained with optimized algorithms [36].

In this research we investigate the potential of CS combined with a low-rank (LR) constraint and sensitivity encoding (SENSE) in achieving fast high-resolution 2D and isotropic 3D 1H-MRSI at 7T for a whole slab of human brain. The ability to accelerate single slice 2D MRSI using CS-SENSE reconstruction was demonstrated at 9.4T by Nassirpour et al. for an FID phase encoded sequence [31] and at 3T by Otazo et al. for a spin-echo PEPSI sequence [27] or by Chatnuntawech et al. for PRESS with random spiral sequence [29]. Additionally, specific reconstruction scheme assuming prior knowledge on MRSI can significantly improve the resulting metabolite distribution [37,38]. Furthermore, the spectral-spatial separability exploited by low-rank modelling has been shown to be very effective to improve SNR and remove nuisance signal, enabling whole brain high-resolution 1H-MRSI at 3T [39-41] with most recent development in modelling including a subspace learned by artificial intelligence [42,43]. To achieve the high robustness that is needed for 3D and 2D MRSI reconstruction we combine the three most powerful methods (CS, LR and SENSE) to successfully deal with the challenges of ultra-high field. We build on our previous formulation of CS-SENSE-LR reconstruction developed at 3T [14,44], and we show that even higher spatial resolutions are reachable for brain metabolic mapping at 7T due to larger SNR and increased spectral dispersion. In addition, we also extend our previous 2D CS-SENSE-LR MRSI to 3D.

Methods

FID-MRSI sequence

Spectroscopic imaging was acquired with a 1H-FID-MRSI [9,10] sequence implemented on a whole-body 7T MRI Magnetom scanner (Siemens, Erlangen, Germany) with a 7T-SC72CD gradient system of 70 mT/m total gradient strength and 200 mT/m/s nominal slew rate, using a 31-channel receive/ birdcage transmit coil and running VB17 software. A slab selective excitation pulse of 1 ms was optimized with a Shinnar-LeRoux algorithm [45] to produce a 6.5 kHz bandwidth and was preceded by four-pulses WET [46] water suppression. The acquisition delay or echo time (TE), between the excitation and the signal acquisition was 1:3 ms in 2D and 0:9 ms in 3D. The free-induction decay (FID) was acquired with 1024 points and 8 kHz sampling rate (spectral window 26.93 ppm), which was followed by spoiler gradients (Fig. 1). The repetition time (TR) was 210 ms and the excitation flip angle was set to 15 degree to prevent T1 weighting in the metabolite signal, considering that the maximum metabolite longitudinal relaxation at 7T in the brain is 1800 ms [47]. K-space was acquired by Cartesian elliptical phase encoding, with both fully sampled and undersampling schemes.

Fig. 1.

(a) Sketch of the FID-MRSI sequence for 2D and 3D. (b) Example of 2D and 3D undersampling patterns in Fourier domain.

Acquisition protocol

An anatomical 3D T1-weighted MEMPRAGE [48] volume was acquired for positioning of the 2D or 3D MRSI. The 2D MRSI was acquired with a 10 mm-thick slice and 210 × 160 mm field of view (FOV). Encoding matrix was 84 × 64 yielding a 62:5 μl voxel volume (2:5 × 2:5 × 10 mm). 3D MRSI was realized with an excited slab of size 210 × 160 × 50 mm (A/P-R/L-H/F). The phase encoding FOV in the head-foot direction was set slightly larger (oversampling) to 60 mm than the 50 mm excited slab to prevent aliasing in this direction due to RF pulse excitation profile. The encoding matrix was set to 64 × 48 × 18 resulting in a 36:5μl voxel volume (3:3 mm isotropic resolution). For signal referencing and to determine the coil sensitivity profiles, unsuppressed water data was acquired with the same FID-MRSI sequence and parameters but with lower resolution: 10 × 10 × 10 mm for 2D and 11:5 × 11:5 × 5 mm for 3D. The FOV size, the excited-slab thickness, the TR and flip angle were identical to main acquisition and were always acquired with a full elliptical encoding.

Undersampled phase encoding

MRSI data were encoded randomly and in a sparse manner over the 3D/2D spatial Fourier domain to grant CS acceleration. Undersampling of the Cartesian encoding was performed by omitting encoding step during the sequential acquisition. The omitting pattern was computed individually in the preparation step of the running sequence. Defining the Fourier domain radius , (k = 0 in 2D), the random sampling was constructed with a density distribution following q−1. The center of the Fourier space with was kept fully sampled (Fig. 1). The size of the fully sampled center was chosen empirically as it represents a reasonable trade-off between reconstruction quality and acquisition time.

MRSI data processing and reconstruction

The removal of lipid signal and residual water was performed as previously described [14]. In short, water signal remaining after the WET suppression was cleared in MRSI raw data of each coil element using the Hankel singular value decomposition (HSVD) method [49]. Afterwards, lipid suppression by metabolite-lipid orthogonality was applied to remove lipid signal in MRSI raw data for each coil element separately [14].

Following the principles of CS acceleration, randomly undersampled MRSI data must be reconstructed with a model that imposes sparsity priors while preserving fidelity with the acquired data [23]. For this study we employed a SENSE reconstruction model for arbitrary trajectory [50,51] combined with total generalized variation (TGV) regularization that imposes sparsity in 1st and 2nd order spatial derivatives. Previous implementation of CS in MRSI also used a combination of Debauchies wavelet and total variation (TV) [28,31] or TV combined with SENSE [29]. In addition, MRSI data were assumed to be low rank and partially separable into spatial and temporal components to enhance SNR [52]. For 2D, we used the model previously published in [14,44] and for 3D, the same model was extended to one extra spatial dimension. We describe this model shortly hereafter. MRSI raw data measured by phased-array coil element c = 1, …, N at time t and at Fourier coordinate k are expressed by the forward model with the integration of the 3D spatial coordinates r over Ω, the object spatial support. The integrand is composed of the transverse magnetization , the coil sensitivity profiles and the spatial frequency shift with ΔB0(r) the field in-homogeneity map in Hz. The transverse magnetization is assumed to be low-rank and separable, i.e. it can be partially separated into K spatial and temporal components, U(r), V(t) [52]:

Combining (1) and (2), the forward model in vectorial notation reads where s is a multi dimensional array containing the 2 or 3 spatial dimensions, the temporal dimension and one coil indexing dimension., and are the Fourier transform operator, the operator applying coil sensitivity profiles and the B0-inhomogeneity correction operator, respectively.

The reconstruction aims to retrieve the spatial and temporal components, U(r), V(t), from the sparsely sampled MRSI data s(k, t). C(r) and ΔB0(r) were computed from water reference acquisitions with the coils sensitivity profiles estimated using ESPIRIT [53] and B0-inhomogeneity profile estimated using multiple signal classification algorithm (MUSIC) [54] on the coil combined water signal spatially interpolated to the metabolite acquisition resolution.

The raw data s cleared of water and lipid contamination (description above) are reconstructed with a low-rank TGV model [44,14]. The spatial and temporal components were then determined by the minimization problem including TGV spatial regularization [55] In the 2D MRSI case, a Hamming filter was integrated in the data fidelity term of the reconstruction model. The modified model eases the convergence of the reconstruction for 2D data with lower SNR while preserving effective voxel size ot the acquistion resolution. A detailed description is given in supplementary material with some comparison on simulated, in vivo and phantom data (Fig. S4). The reconstruction was performed in Matlab (The Math-Works, Inc., Natick, Massachusetts, US), required a minimum of 64 RAM GB and necessitate 12hour computation time for the 3D MRSI using 8 cores of a 3.00 GHz Intel(R) Xeon(R)-E5 CPU.

LCModel quantification

Metabolite concentrations were quantified using LCModel [56] to fit the MRSI dataset = UV that resulted from the optimization process in (4). A reference basis was simulated using GAMMA package [57] with acquisition parameters identical to the MRSI sequence.

The water reference data underwent the same reconstruction as the MRSI dataset and the reconstructed water signal was used as reference by LCModel. Thanks to the ultra-short TE and the low flip-angle, no T1 or T2 relaxation correction are necessary for metabolite quantification relative to the internal water reference signal.

The following metabolites were included in the simulated basis: N-acetylaspartate (NAA), N-acetyl aspartylglutamate (NAAG), creatine (Cr), phosphocreatine (PCr), phosphorylcholine (PCh), glycerophosphorylcholine (GPC), myo-inositol (Ins), scyllo-inositol, glutamate (Glu), glutamine, lactate, beta-glucose, alanine, taurine, aspartate, gamma-aminobutyric acid and glutathione. The removal of lipid signal by orthogonality may result in a strong baseline distortion under NAA singlet peak at 2 ppm in the form of a negative undulation. To help the LCModel quantification to cope with this distortion, a 20 Hz-broad peak at 2 ppm and with negative phase was added to the basis. LCModel control parameters are given in the supplementary material.

The results of LCModel fitting for each voxel were further used to create spatial maps for the concentration of each metabolite. The spectral quality was assessed through the goodness of the fit and was estimated by the residuals root mean square (RMS) for each voxel.

High-resolution structural and metabolic phantom

To assess precision and accuracy of the CS-SENSE-LR FID-MRSI, a high-resolution structural and metabolic phantom containing tubes of several diameters was measured with 2D FID-MRSI without CS acceleration and with a 2 mm in plane resolution. To simulate acceleration, the data were undersampled retrospectively with acceleration factor 2,3,4 and 5 following the same probability distribution in the Fourier domain as the accelerated acquisition. For comparison, a reconstruction by the forward-model adjoint operator (Fourier transform and coil combination, details in supplementary material, Model Alternatives) was also performed on fully sampled data. The geometry and molecular contrast of our custom made phantom is similar to Derenzo phantom [58,59] used for quality control and quantification in PET molecular imaging. Our custom made phantom consists of a large cylindrical container of clear cast acrylic material (outer diameter (OD)= 15.24 cm, inner diameter (ID)= 13.33 cm, Mc-Master-Carr), 5 sets of tubes corresponding to the diameters 2, 4, 6, 8 and 10 mm, and a tube holder of size equal to the ID that is firmly fixed on the inner wall of the container. For each diameter size the set consisted of 6 tubes arranged. The tubes fixed by the holder are separated by a distance equal to twice the inner diameter of the tubes in a triangular close packed configuration as shown in Fig. 2. The tubes were filled with specific metabolite solutions of six different concentrations and mixtures based on their position in each size set (1–6) (top right table Fig. 2). Magnevist (Gd-DTPA) was added (1 mL/L) in each tube to shorten T1 and create T1-weighted contrast for structural MRI. The whole tube structure was inserted in the large container which was filled with 10 mM NaCl solution. The tubes and holder of the phantom were designed using 3D computer-aided design (CAD) software (rhinoceros 6.0, Robert McNeel & Associates) and fabricated using 3D printer (Formlabs, Form 2, Somerville, MA USA). The printing material was clear resin (RS-F2-GPLC-04) which is a mixture of acrylated monomers, acrylated oligomers, and photoinitiators.

Fig. 2.

2D metabolite maps of the high structural and metabolic phantom acquired with 2 × 2 × 10mm3 resolution and accelerated retrospectively with factor 1, 3 or 5 (left). The actual concentration in each tube is given in the top table in mM. The tube labeling is the same for all five sizes. Right bottom, voxel-average metabolite ratio over creatine in tubes 1 to 6 of all sizes as function of the actual concentration ratio at different acceleration factors. The error bars represent the standard deviation over all the voxels within a certain tube size.

Healthy volunteer MRSI measurement

To illustrate the feasibility of the CS-SENSE-LR FID-MRSI, an acquisition protocol with 2D MRSI acceleration factor 2:5 (6 min) and 3D MRSI acceleration factor 4:5 (20 min) was acquired over five young healthy volunteers. In addition and to measure the effect of the CS acceleration by retrospective undersampling and reproducibility regarding the accelerated acquisition, a highly sampled 3D MRSI with acceleration factor 2:5 (34 min) was acquired in volunteer 2 and a fully sampled 2D MRSI (15 min) was acquired for three volunteers. The protocol was approved by the institutional ethics committee and written informed consent was given by all subjects before participation.

Results

High-resolution phantom

The metabolite images of the high-resolution phantom reconstructed with the CS-SENSE-LR model display the tube cross-sections that are distinguishable for all diameters in comparison to the adjoint operator reconstruction that does not always permit the distiction of the 2 mm tubes (Fig. 2). Patterns of concentrations corresponding to the metabolite content are visually observable for NAA and Cho in tubes with diameter 10, 8 and 6 mm. The Cre images and the signal variability among tubes containing all the same concentration indicate the difficulty to obtain homogeneous signal throughout the phantom for the adjoint or the CS-SENSE-LR reconstructions. This is probably caused by the strong B0 inhomogeneities and the difficult correction by water reference signal with the presence of the resin (tube construction material) that creates signal drop around the tubes. As a consequence, the water signal could not be used as reference for absolute quantification but instead the ratio to Cre were studied. In Fig. 2 right, the measured concentration ratios versus the actual concentration ratios are displayed. The measured concentration values correspond to the voxel average in tubes from size 2 to 10 mm at one position (1–6, as label in top of Fig. 2) with the standard deviation as error bar. Accuracy of the results can be assessed by the distance to the diagonal (dash blue line) and the correlation coefficient square R2. NAA/Cre, Cho/Cre and Glu/Cre show good agreement between measured ratio and actual ratios for both the CS-SENSE-LR and the adjoint reconstructions although low concentration of Cho or Glu seems to be overestimated. The results are practically unchanged when using the CS-SENSE-LR model and increasing acceleration factors, which supports the accuracy of the reconstruction. Further comparison between the Adjoint, the CS-SENSE-LR and other model alternative are presented in supplementary material and Fig. S3.

Retrospective acceleration on healthy volunteer MRSI measurement

To evaluate the performance of CS acceleration on reconstructed metabolite maps and volumes in volunteers, 3D FID-MRSI data were reconstructed with several acceleration factors resulting of retrospective undersampling. The effect of CS acceleration are illustrated on 3D metabolite volumes of one volunteer in Fig. 3 in an orthogonal view. While acceleration of the MRSI data was progressively increased up to 6.5, only a slight loss of fine detailed contrast can be observed. A normalized root mean square error (NRMSE) over the brain for each metabolite is computed for each acceleration factor relative to the no-acceleration dataset. Let being the metabolite concentration at voxel υ for CS acceleration factor x, the NRMSE for the metabolite map reads . The NRMSE values show a linear behavior with the acceleration factor in agreement with previously published work [60,61,31]. The spectral quality illustrated by the three sample spectra remains high for all acceleration. This is supported by the constant fitting residuals RMS as function of the acceleration factor shown in Fig. 3. This results are consistent with similar retrospective acceleration analyses performed on 2D FID-MRSI and simulation data in supplementary material Fig. S1 and S7.

Fig. 3.

Top, 3D 3:3mm isotropic metabolite volumes acquired with FID-MRSI and different acceleration factors. Bottom left, T1-weighted images of the corresponding metabolite slices in the three direction are shown next to three sample spectra at acceleration factor 2:5, 4:5 and 6:5. Bottom right, as indicator of spectral quality, the residuals root mean square (RMS) of the LCModel fit is displayed versus the acceleration. Also, the normalized root mean square error (NRMSE) of each metabolite over the whole slab relative to A.F.=2.5 versus the higher acceleration factors. Similar retrospective acceleration analysis for 2D dataset can be found in supplementary material Fig. S1. Metabolite decomposition of the voxel 1 fitting is presented in supplementary material Fig. S6.

The 2D metabolite maps reconstructed from the five volunteer measurements with acceleration factor 2.5 are shown in Fig. 5 for NAA, tCr, Cho, Glu, Ins and NAAG with the T1-weighted image corresponding to the slice location. The maps reveal spatial anatomical patterns that are specific for each metabolite distribution in the brain and that are common to all five volunteer datasets. NAA distribution is homogeneous throughout the slice whereas tCr concentration is higher in grey matter (GM) than white matter (WM). Cho map exhibits a major concentration in the frontal WM while concentrations are low in the occipital lobe. The strongest GM/WM contrast is present in Glu concentration maps where the cortex can be distinguished. NAAG is clearly only detected in WM although it is the least accurate map due to NAAG’s low signal and overlap with NAA. The same spatial features characterizing metabolite distributions can be observed in the 3D metabolite volumes. These were measured on the same volunteers with 3D CS-SENSE-LR FID-MRSI accelerated with acceleration factor 4:5 and one volunteer’s data are presented in Fig. 4. Results of two additional volunteers are shown in the supplementary material Fig. S2. The T1-weighted images show the anatomical content of the slab. A quantitative atlas-based analysis was performed to estimate the metabolite concentrations in the GM and WM of each cerebral lobe (detailed description of the model in supplementary material, Atlas-Based Analysis). The WM to GM ratio of concentrations was computed in all atlas structure and reported in Table 1 for the 2D and 3D healthy volunteer data. The tCr and Glu ratios exhibit values lower than 1 that illustrate a constant high concentration in GM than WM, in agreement with the observation described above. In contrast, Cho ratio values are greater than 1 and NAA values are approximately equal to 1. The NAAG values are particularly high because this metabolite is mainly located in WM, in agreement with previous literature results [62,63]. Metabolite ratios to tCr were also computed using the atlas-based analysis. The results presented in Table 2 are in good agreement with previously reported values measured at ultra-high field with FID-MRSI [10,64]. The strong B1 and B0 inhomogeneities resulting in intensity variation and frequency shift were corrected by the model reconstruction using the water-unsuppressed measurement (supplementary material Fig. S10). An estimate of the reproducibility and effect of acceleration was performed by comparison between measurements of the highly sampled data and accelerated acquisition in 2D and 3D. The results that are presented in supplementary material Fig. S8 and S9 show good agreement between the two consecutive measurements and the error maps remain under 25% value. The RMSE in 2D and 3D compare to the values reported for the retroscpective acceleration in Fig. 3 and Fig. S1, and suggest that the acceleration is the major cause of error. The systematic higher RMSE values for volunteer 4 compared to volunteer 1 tend to indicate slight motion between the measurement.

Fig. 5.

2D CS-SENSE-LR FID-MRSI metabolite maps of five healthy volunteers measured with acceleration factor of 2:5 for a total acquisition time of 6 min. The colour scale goes from 0 till the 95th percentile for each metabolite separately. The T1-weighted anatomical image corresponding to the MRSI slice is shown to the right.

Fig. 4.

3D CS-SENSE-LR FID-MRSI metabolite volumes of one healthy volunteer performed with acceleration factor of 4:5 for a total acquisition time of 20 min. The colour scale goes from 0 till the 95th percentile for each metabolite separately. The T1-weighted anatomical images to the right show the MRSI slab positions. 3D CS-SENSE-LR FID-MRSI metabolite volumes of two additional volunteers can be found in the supplementary material Fig. S2.

Table 1

Ratios between the WM and the GM concentrations in each anatomical lobe and average accross all volunteers. In 2D values for the three lobes that were consistently covered are reported.

WM to GM ratio in anatomical lobes
3D	NAA	Cr + PCr	Ins	PCh + GPC	Glu	NAAG
Frontal	0.99	0.89	0.99	1.08	0.82	1.76
Limbic	0.92	0.81	0.93	1.00	0.73	1.81
Parietal	0.92	0.85	1.01	1.08	0.78	1.95
Occipital	1.02	1.01	1.17	1.17	0.90	1.47
Temporal	0.86	0.81	0.994	1.03	0.75	1.47
2D	NAA	Cr + PCr	Ins	PCh + GPC	Glu	NAAG
Frontal	1.02	0.94	1.14	1.17	0.82	2.91
Limbic	0.99	0.88	1.13	1.12	0.78	2.82
Parietal	0.95	0.86	1.03	1.09	0.74	1.62

Table 2

Metabolite concentrations relative to Cr + PCr in each anatomical lobe and average accross all 3D volunteer measurements.

Metabolite ratio to Cr + PCr in anatomical lobes and tissues
3D	NAAtCr		Ins tCr		PCh+GPCtCr		Glu  tCr		NAAtCr
	GM	WM	GM	WM	GM	WM	GM	WM	GM	WM
Frontal	1.75	1.93	0.81	0.91	0.29	0.35	1.08	0.99	0.17	0.34
Limbic	1.67	1.90	0.82	0.95	0.30	0.37	1.10	0.98	0.15	0.35
Parietal	1.69	1.85	0.78	0.93	0.26	0.33	1.13	1.03	0.13	0.32
Occipital	1.65	1.68	0.75	0.88	0.26	0.30	1.09	0.98	0.18	0.26
Temporal	1.65	1.75	0.81	0.99	0.30	0.37	1.05	0.97	0.16	0.31

Discussion

Brain high-resolution 2D/3D MRSI at ultra-high field can be performed with greatly reduced acquisition time and high sensitivity by combining FID-MRSI acquisition with CS-SENSE-LR acceleration and reconstruction. This achievement was enabled by the use of several technical improvements. The TR of the FID-MRSI sequence was minimized to 210 ms. Cartesian phase encoding was accelerated by random undersampling. The CS-SENSE-LR model enables accurate MRSI data reconstruction while enhancing signal with LR constraint.

The results of retrospective acceleration applied on the FID-MRSI data in 3D and 2D (Fig. 2, Fig. 3 and supplementary material Fig. S1 and S2) showed that the CS-SENSE-LR model allows for accurate reconstruction albeit random undersampling is performed by factor 3 or more. For 3D MRSI (Fig. 3), the reference dataset was not fully sampled but was acquired with a 2.5 acceleration factor. A fully sampled 3D acquisition would be lengthy (85 min) with a high risk of head motion and large scanner frequency drift. It was assumed that considering the CS results in 2D (Fig. 2, supplementary Fig. S1 and [14]), an acceleration factor of 2.5 in 3D should contain minimal acceleration distortion and represent a good starting point for retrospective acceleration.

Optimal accelerated acquisition

Acceleration factor of 3 in 2D (5 min) and 4:5 in 3D (20 min) were found to be the best compromise with a strong reduction in acquisition time while presenting minimum effects on metabolite maps. The optimal acceleration factor of 3 found here for 2D and 2:5 mm in-place resolution compare to the factor of 4 in 2D and 3:1 mm resolution at 9.4T found in [31]. The difference might be explained by higher SNR at 9.4T but also by the qualitative determination of the factor. 2D MRSI results may suffer from lower SNR in comparison to 3D MRSI even though voxel volume is larger in 2D case. This is a direct consequence of the larger volume excited by the 3D sequence and might favor a 3D over 2D multislices MRSI acquisition for particular application that requires high SNR for detection of low signal metabolites. Compared to our previous work [14], we were able to increase the matrix size and reduce acquisition time in 2D MRSI from an resolution of 3:3 × 3:3 × 10 mm3 in 11 min at 3T to 2:5 × 2:5 × 10 mm3 in 6 min at 7T thanks to higher SNR and shorter possible TR. The accelerations obtained with CS for our 2D protocol are similar to accelerations obtained for similar protocols by parallel imaging [19,65], with the added benefit that CS-SENSE-LR reconstruction provides effective spectral denoising and less structured artifacts from undersampling. In addition, the reconstruction presented here has the advantage to produce metabolite image with an effective voxel size identical to the nominal size for fully sampled dataset without acceleration as illustrated with the high resolution structural phantom (Fig. 2) and the simulated data (Fig. S4). When acceleration by random undersampling is performed, a loss of resolution can be visible and compare to the g-factor noise increase for GRAPPA acceleration. The reproducibility results show good agreement between the two consecutive acquisitions and consistent results in 2D and 3D, but with slightly higher RMSE in 2D (supplementary material Fig. S8 and S9). The lower SNR in 2D acquisitions might explain partially the lower reproducibility.

The original 2D method was also successfully extended to 3D here while benefiting from 7T advantages to reach an acquisition resolution of 3:3 mm3 in 20 min. Our 3D protocol using CS to acquire a matrix of 64 × 48 × 18 in 20 min compares well with a recently published 3D protocol [64] using concentric rings to acquire a matrix of 80 × 80 × 47 in 15 min or 64 × 64 × 39 in 9 min.

The TGV regularization parameter was adjusted to λ = 1 × 10−3 for the 2D reconstruction (same value as for 3T study [14]) and λ = 3 × 10−4 for the 3D reconstruction. These same values were observed to be optimal for all volunteers in 3D or 2D as shown in supplementary material (Fig. S5). Therefore, no adjustment of the regularization parameter is necessary for each subject but should be slightly readjusted if acquisition parameters such as slab thickness, resolution, flip-angle or coil setup are modified. However, this is generally true not only for our reconstruction but for all model based reconstructions, including LCModel fitting (e.g. a change in sequence or B0 field requires a change in fitting basis set). The number of components in the reconstruction (K in (2)) was chosen as follow. As described in [14], the initial estimate of the spatial and temporal components are computed by SVD on the adjoint solution. Initial spatial and temporal components are then reviewed and K was qualitatively chosen as being the minimum number of components containing some signal distinguishable from noise. K was set to 26 for 2D and 40 for 3D.

Artifacts and distortions

For the highest acceleration (≥ 4 in 2D and ≥ 6 in 3D), a loss in fine details of the metabolite maps is observed. The signal containing these sharp contrasts is stored in the high spatial frequencies located in the outer Fourier space, precisely where most of the sampling points are removed. It is therefore consistent that with utmost undersampling, small contrast features are lost. This finding is consistent with the literature [51,25] where strong CS MRI acceleration results in noticeable blurriness or loss of image resolution. These image alteration at high CS acceleration can affect small features or large structure differently and therefore cannot be described as a unique increase of the effective voxel size.

Nassirpour et al. observed an increasing contamination by lipid signal for higher acceleration [31] but that was not the case in our results and is probably due to the fact that in our proposed pipeline, lipid and water suppression take place before the reconstruction whereas in [31] these are performed after reconstructing the MRSI. In particular, in the previous work [31] which investigated separately the performance of low-rank and CS reconstruction, it was observed that low-rank method provided less lipid artifacts while CS provided lower CRLB for metabolite fitting. Hence, by combining LR and CS in our reconstruction we can take advantage of both properties to improve data quality and quantification. Spectral quality observable in Fig. 3 with the different spectra and assessed by the fitting residuals is unaffected by CS acceleration, in agreement with previous observations [65].

strong inhomogeneities present at ultra-high field were implicitly corrected by referencing metabolite signal with unsuppressed water data. Frequency shift yielded by B0 variation were corrected within the CS-SENSE-LR model based on a B0 fieldmap computed from the unsuppressed water data but apparent loss of signal may still occur in region where local B0 inhomogeneity is large and reduces strongly the metabolite .

The lipid suppression by metabolite-lipid orthogonality permits a complete removal of lipid contamination from the 2D and 3D FID-MRSI datasets. However, a consecutive distortion of the spectra baseline under the NAA singlet at 2 ppm may be observed. The sample spectra shown in Fig. 3 display this distortion with a broad baseline dip resulting in a ‘W’ shape at 2 ppm. This baseline deformation was circumvent with the introduction on extra peak in the LCModel basis at 2 ppm. This peak was introduced with a negative phase opposed to NAA singlet and with a broad 20 Hz width to prevent cross-correlation with NAA during fitting by LCModel. As results, LCModel succeeds in close fitting of the data with a smooth baseline as shown in Fig. 3 and supplementary material Fig. S6. This effect of the lipid suppression on spectra is an important point. Although the additional peak added to the basis copes correctly with the baseline pit, we cannot fully discard an effect on NAA quantification. This point shall be addressed more in details in a future publication. The computation time necessary for the reconstruction of 3D high-resolution FID-MRSI data that was about 12 h and might be a limitating factor in the application of the technique. However, we expect that with a reconstruction code fully optimized and with specific dedicated computational resources such as graphics processing units, it could be computed in significantly shorter time such in the case of MRI [36].

Literature background

To the best of our knowledge, the first 3D images of metabolite distributions over extensive brain coverage were published by Duijn et al. [66] with low-resolution Cartesian encoding at 2T MRI scanner. These results were followed by development of fast encoding EPSI/PEPSI techniques [67-70], spiral encoding [71,72] and rosette spectroscopic imaging [73]. Group analysis of whole brain 3D MRSI provided evidence of tissue and metabolic specific differences across the brain [74,75]. Recent publication of metabolite distributions at ultra-high field using concentric-ring trajectories showed an apparent GM/WM contrast in high resolution for Cho/tCr and Glx/tCr volumes [64,21]. To highlight the sensitivity and the reproducibility of the method presented here, these metabolite features were retrieved in both 2D and 3D maps (Fig. 4 and 5, Table 1) with higher contrast of GM compared to WM in tCr and Glu maps. The reconstructed NAAG maps show a typical presence in central WM as previously shown in literature [62,63]. The Cho distribution exhibits a characteristic presence in WM with high signal in frontal lobe and low level in occipital lobe.

Additional acceleration of our 3D CS-SENSE-LR MRSI method could be achieved by combining CS with non-cartesian trajectories as shown at 3T for MRSI with PRESS excitation [29]. Also, the combination of non-cartesian encoding with low-dimensional models was demonstrated to be particularly efficient with SPICE [39,43] or LR reconstruction of spiral trajectories [40].

Conclusion

An acquisition and reconstruction combination was proposed for CS accelerated high-resolution 1H-FID-MRSI at ultra-highfield. The reconstruction model that includes LR and TGV constraints, enables acceleration by phase-encoding random undersampling which is demonstrated for the first time in 3D over human brain, with little effects on metabolite maps. This MRSI method grants measurement of high-resolution whole-slab metabolite distributions in short acquisition times: 5 min for 2D MRSI 2.5 mm resolution and of 20 min for 3D MRSI at 3.3 mm isotropic.