Retrospective frequency drift correction of rosette MRSI data using spectral registration

Frequency drift correction is an important postprocessing step in MRS that yields improvements in spectral quality and metabolite quantification. Although routinely applied in single‐voxel MRS, drift correction is much more challenging in MRSI due to the presence of phase‐encoding gradients. Thus, separately acquired navigator scans are normally required for drift estimation. In this work, we demonstrate the use of self‐navigating rosette MRSI trajectories combined with time‐domain spectral registration to enable retrospective frequency drift corrections without the need for separately acquired navigator echoes.


INTRODUCTION
Proton ( 1 H) MRSI enables mapping of tissue metabolite concentrations in two or three dimensions in vivo. [1][2][3] Over the past several decades, development of MRSI acquisition, reconstruction, and data-processing techniques has led to improved spatial coverage and resolution, acquisition speed, artifact rejection, and metabolite quantification. [4][5][6][7] In particular, the development of spatial-spectral encoding techniques such as echo planar spectroscopic imaging, 8,9 spiral, 10 concentric rings, 11 radial, 12 rosette, 13 and stochastic trajectories 14 has enabled rapid and high-resolution MRSI data acquisition.
In vivo MRS/MRSI experiments are negatively affected by temporal fluctuations in the main magnetic field (B 0 ). Drifts in the B 0 field primarily occur as a result of heating of passive shim elements by the gradient coils, and subject motion (both physiological and bulk). 15 This results in shot-to-shot accrual of frequency and phase variations between repetitions, thereby causing incoherent averaging and/or encoding, which in turn results in spectral peak broadening, lineshape distortion, spectral misalignment, and signal-to-noise ratio (SNR) losses in the resulting spectra.
In single-voxel MRS, B 0 drift can be corrected retrospectively by aligning each acquired average to a reference scan (usually one of the averages in the series). Several alignment approaches have been proposed, including correlation-based methods, 16 residual water peak correction, 17 metabolite lineshape fitting, 18 and time-domain spectral registration (SR). 19 However, drift correction in MRSI proves to be more challenging because of the sampling of different k-space coordinates across repetitions. If the resulting repetitions do not share any common k-space points (as is the case in most MRSI sampling trajectories), then spectral alignment cannot be used for drift estimation. In such cases, techniques using separately acquired navigator echoes [20][21][22][23] have been proposed for drift correction in MRSI. However, this approach complicates implementation, and may also result in longer acquisition times.
In contrast, centric k-space trajectories such as spiral, radial, or rosette trajectories do contain a common k-space point across all repetitions, namely, the center of k-space or k = 0. Such centric trajectories are said to be "self-navigating." In this paper, we demonstrate the use of the self-navigating rosette MRSI trajectory for characterizing and correcting for scanner B 0 drift using time-domain spectral registration in vivo. Specifically, from each repetition (i.e., each TR cycle), we isolate the FID from the center of k-space (i.e., the k = 0 FID). Then, using the SR algorithm, the k = 0 FIDs are retrospectively aligned to a reference scan (i.e., the first k = 0 FID in the series) to estimate the frequency offset at each repetition. The estimated frequency offsets are then used to apply corrections throughout k-space.
To demonstrate the efficacy of this approach, rosette MRSI data were acquired in 5 healthy volunteers at 3 T.
To assess spectral quality, the linewidth (LW) and SNR of the N-acetylaspartate (NAA) peak were assessed throughout the brain before and after drift correction. Next, we assessed the improvement in quantification precision of brain metabolites obtained using LCModel (Stephen Provencher, Inc.). 24 Results from this study demonstrated a significant improvement in spectral quality following drift correction.

Study participants
Rosette MRSI data were acquired in 5 healthy volunteers (4 males, 1 female; age: 29 ± 6 years) with no history of brain disease. All study procedures were approved by the Research Ethics Board of Sunnybrook Hospital.

Data acquisition
A 2D PRESS-localized rosette MRSI sequence ( Figure 1A) was implemented on a 3T Prisma MR scanner (Siemens) with a 20-channel phased array head coil and the following scan parameters: TR/TE = 1200/30 ms, flip angle = 72 • , 8 averages, PRESS excitation region = 90 × 100 × 15 mm 3 , encoding FOV = 240 × 240 mm, matrix size = 48 × 48, nominal voxel size = 5 × 5 × 15 mm 3 , spectral bandwidth = 1587 Hz, and slice-selective RF pulse center frequency = 3.0 ppm. The rosette trajectory ( Figure 1B) was implemented with 126 k-space samples per rosette petal, 576 FID points per k-space position (i.e., each rosette petal was traversed 576 times sequentially), and 76 rosette petals/shots. During the acquisition of each petal/shot, the ADC was open for a total duration of 362.88 ms, and a total of 72 576 sample points were acquired. Adjustment of all first-order and second-order shims was performed using the Siemens gradient-echo shim method (i.e., the "Brain" shim mode). 25 Before excitation, water suppression enhanced through T 1 effects 26 was performed, consisting of three frequency selective pulses each with a duration of 25.6 ms and a bandwidth of 80 Hz. Four outer-volume suppression bands of 30-mm thickness each were manually placed around the edges of the PRESS localization region to further suppress extracranial lipid signals. The sequence was implemented in long-term averaging mode, whereby all 76 rosette petals are traversed before repeating signal averages. The total scan time was approximately 12 min. Before the rosette MRSI scan, a 3D T 1 -weighted, MP-RAGE image (1-mm isotropic resolution) was acquired for localization purposes. 27 The MRSI slice was oriented axially and positioned just above the lateral ventricles ( Figure 2).
Prior to implementing rosette MRSI in vivo, phantom experiments (using a GE Braino phantom; GE Healthcare) were performed using the same scan protocol as above, for validation purposes ( Figure 1D).
The MAGNETOM Prisma scanner can achieve a maximum gradient amplitude of 80 mT/m and a slew rate of 200 T/m/s. Using a spectral bandwidth ( ) of 1587 Hz yielded a peak gradient strength of 11.71 mT/m along the x-direction and y-direction, and 7.28 mT/m along the z-direction. The peak slew rate in the x-direction and y-direction was 127.97, and 36.66 T/m/s along the z-direction. Data were extracted from the scanner console in TWIX (.dat) format, which contains individual raw FIDs for each RF channel and repetition.

Data processing
Data processing and reconstruction were performed in MATLAB (MathWorks, Natick, MA, USA) using the FID-A Toolkit. 28 First, k-space density compensation was performed using the Voronoi density estimation method. 29,30 A nonuniform discrete Fourier transform was then performed along the spatial dimensions, following by RF channels being combined as follows: At each voxel position, the coil's phases were first adjusted to match the phase of the first coil; then, sensitivity weights for coils were estimated by the formula where S is the intensity of the first point in the FIDs, and i is the coil number. Coils were then summed by the equation F = sum(W(i) * f (i)), where W(i) is the ith coil's sensitivity weight; f (i) is the ith coil FIDs; and F is the final summed signal. 31 This way, signals from multiple coils were combined for optimal SNR at every image voxel. The phase and weight maps calculated for the coil combination of water unsuppressed data were applied to the water-suppressed data, after which spectral Fourier transform was performed to obtain the resulting 48 × 48 MRSI grid. Residual water peak removal was implemented using L2 regularization. 32

Frequency and phase drift correction
Frequency drift correction was performed prior to the density compensation in the processing pipeline described above. k = 0 FIDs were isolated from each repetition by extracting the first and each subsequent 126th point of each ADC readout (126 being the number of ADC points per k-space loop) ( Figure 1C). Using the first k = 0 FID as a reference, frequency and phase offsets for each subsequent repetition were estimated using time-domain spectral registration, which uses nonlinear least-squares minimization to estimate the frequency and phase shift required to optimally align each FID signal to the reference spec- This minimization was implemented in MATLAB using the "nlinfit" function provided by the MAT-LAB Statistics toolbox, which returns the least square parameter estimates based on the Gauss-Newton algorithm with Levenberg-Marquardt modifications for global convergence. Because the early part of the FID always contains the highest signal amplitude, the SR algorithm uses only the early part of the FID for this alignment procedure. Specifically, the reference FID (R(k = 0, t)) and alignment FIDs (S n (k = 0, t)) were truncated at t = 500 ms. Once the frequency and phase offsets were estimated for each k = 0 FID, the corresponding corrections were applied to the full ADC signals to achieve a frequency and phase drift correction across the whole of k-space (not only the k = 0 extracted FIDs). This was achieved by multiplying the FIDs acquired at all k-space positions along each nth rosette petal, S n (k, t), by a frequency and phase correction term as follows: where f opt n and opt n are the optimized frequency and phase correction terms, respectively, determined by alignment of the nth k = 0 FID to the reference k = 0 FID. Estimation and correction of drifts were performed individually for each RF channel. The MATLAB code used for spectral registration is provided in the FID-A repository on GitHub (op_CSISpecReg.m). 28 The performance of time-domain spectral registration was assessed by visual inspection of frequency and phase drift curves, and by assessment of spectral quality (SNR and LW) before and after correction. SNR was defined as the amplitude of the largest metabolite peak (NAA) in the absolute part of the spectrum divided by the SD of the real part of the spectrum in a region with no spectral peaks (−2.0 to 0 ppm), and spectral LW was defined as the FWHM of the NAA peak in the water-suppressed spectra. SNR and LW maps were generated in FID-A. Comparisons of SNR and LW before and after field drift correction were performed using a two-tailed Student's paired t-test, with a probability of p < 0.05 considered as statistically significant.

Metabolite quantification
Processed MRSI data were fit using LCModel. 24 A non-water-suppressed rosette MRSI acquisition (two averages) was used to provide water referencing on a voxel-wise basis in LCModel. The LCModel default values for WCONC, ATTH2O, and ATTMET were used, and no additional correction for T 1 -relaxation or T 2 -relaxation of water or metabolites was performed.
The non-water-suppressed scan was also used for eddy current correction. An LCModel basis set was simulated in FID-A and included 21 brain metabolites: alanine  .75 ppm). The parameters (frequency and LW) of these MM components were estimated based on acquired metabolite-nulled localized brain MRS spectra from healthy volunteers at 3 T. Spectra were fitted within the chemical shift range 0.5-4.2 ppm. All spectra from within the excited PRESS region (90 × 100 mm 2 in-plane) were included for further analysis in all subjects, and the following metabolites of interest were considered: tCh (GPC + PCh), tNAA (NAA + NAAG), tIns (Ins + Gly), tCr (Cr + PCr) and Glx (Glu + Gln), as these metabolites consistently showed reasonable precision in their quantification (i.e., Cramér-Rao lower bound [CRLB] < 20%). CRLB values for LCModel quantified data were compared using a two-tailed Student's paired t-test to assess changes in fit uncertainty following drift correction. The details of the rosette MRSI acquisition and analysis are also included in the MRS reporting checklist (Table S1). 34

Field drift characterization
The measured frequency drift from one of the in vivo rosette MRSI acquisitions is shown in Figure 3A. The resonant frequency drifted approximately linearly by a total of 5.7 Hz across all scans during the 12-min acquisition F I G U R E 3 (A) Plot of frequency drift measured for each coil channel and shot over 608 scans (76 shots × 8 averages). The frequency drift for this scan was estimated to be approximately 0.48 Hz/min, with a total drift of 5.7 Hz. (B) Water-suppressed spectra from an arbitrarily chosen gray-matter (top) and white-matter (bottom) voxel in the reconstructed MRSI grid, before (blue line) and after (red line) drift correction. Metabolite peaks appear to have higher amplitude following drift correction.
for the representative data set. The average total frequency drift measured across all data sets was 5.8 Hz (0.48 Hz/min). Because none of the scans in this study were preceded by gradient-heavy EPI or diffusion imaging sequences, 35 the observed frequency drift is likely due to gradient heating induced by the rosette gradients themselves. Figure 3B shows water-suppressed spectra before and after field drift correction from one arbitrarily chosen voxel in gray matter and one in white matter from the reconstructed 1 H-MRSI grid. Visually, the metabolite peaks appear to have higher amplitude following drift correction.

3.2
Assessment of spectral improvement following drift correction Figure 4 shows NAA-SNR and NAA-LW maps measured before and after drift correction in a representative subject. Across all voxels in all subjects, the average values of SNR increased by 12.9% from 29.9 ± 2.6 to 33.8 ± 3.0 [t(528) = −16.05, p < 0.05], and the NAA LWs decreased by 18.5% from 9.7 ± 0.6 Hz to 7.9 ± 0.8 Hz [t(530) = 29.99, p < 0.05]. Figure 5 summarizes the improvement in SNR and LW following drift correction for each individual subject as well as for the entire study cohort. To assess the potential change in metabolite concentration estimates and their uncertainties after spectral registration, metabolite concentration maps were generated for tNAA, tCh, tCr, and tIns based on the LCModel fitting results ( Figure 6). Table 1 illustrates the mean LCModel absolute concentration estimates and the corresponding CRLB values before and after drift correction for the selected metabolites in all subjects as well as the whole study cohort. The mean CRLB values for tCh, tNAA, tIns, and tCr showed a statistically significant decrease (p < 0.05), with reductions of 3.9%, 4.0%, 12.6%, and 5.3%, respectively. In contrast, the mean CRLB value for Glx showed a non-significant decrease of 0.8% (p > 0.05). This may be attributed to the relatively low SNR of the Glx peaks, which made it difficult to detect any statistically significant changes in Glx peak heights or fitting quality after drift correction.

DISCUSSION
This study demonstrates the use of the self-navigating rosette MRSI trajectories to retrospectively correct frequency and phase drifts in the acquired in vivo MRSI data. Drift correction was possible without the need for separately acquired navigator echoes, as the rosette trajectory collects an FID at the center of k-space on every shot, which can be used for frequency alignment. Other centric k-space trajectories such as spiral or radial MRSI would be similarly amenable to this approach. However, approaches using other commonly used MRSI trajectories, including conventional MRSI, concentric rings MRSI, or echo planar spectroscopic imaging, cannot benefit from this approach, as these acquisitions do not share any common k-space points from shot to shot.

F I G U R E 4 (A,B) Illustration of N-acetylaspartate (NAA) SNR (A) and
NAA linewidth (in Hz) (B) maps generated using FID-A, before (left) and after (right) field drift correction in a representative subject. NAA SNR and NAA LW improved significantly (p < 0.05). The improvements were reasonably uniform across the region of interest.

F I G U R E 5
(A,B) Box plot representation of SNR (A) and (B) linewidth (in Hz) (B) for each individual subject as well as for the whole study cohort, before (blue) and after (red) drift correction. The overall SNR increased by 12.9%, and the overall spectral linewidth decreased by 18.5%.
A recent multicenter study investigated the magnitude of frequency drift during single-voxel MRS experiments at 3 T. 36 Across 95 3T scanners from three major vendors (GE, Philips, and Siemens), the average rates of frequency drift were 0.29 and 0.43 Hz/min before and following a gradient-intensive echo-planar

F I G U R E 6
LCModel metabolite maps for the absolute concentrations of tCh (GPC + PCh), tNAA (NAA + NAAG), tIns (Ins + Gly), tCr (Cr + PCr), and Glx (Glu + Gln) in native 48 × 48 resolution before (above) and after (below) field drift correction, shown in a representative subject. The estimated concentrations of tCh, tNAA, tIns, and tCr increased slightly following drift correction, while the overall concentration of Glx remained unchanged.
functional MRI acquisition, respectively. In the current study, although the rosette acquisition did not follow any gradient-intensive imaging sequences, we observed a relatively high drift rate of 0.48 Hz/min (average across 5 participants). This high rate of drift suggests that spatial-spectral MRSI sequences are likely associated with more gradient heating and greater frequency drift compared with single-voxel MRS. This result is not unexpected, given the intensive sinusoidal gradient trajectories used, and suggests that as rapid spatial-spectral MRSI approaches emerge into mainstream use, frequency drift correction will be an important consideration. (Note that in the current study, scans were not performed early in the morning, and thus were likely not affected by drift effects that are known to occur when the system is waking from overnight energy economy mode). The drift rate results obtained in this study were notably lower than the typical drift rates observed in diffusion MRI studies (|f | = 3.1 Hz/min), 37 indicating the relatively lower impact of rosette MRSI on gradient coil heating. Future work should aim to characterize typical field drift across multiple scanners and its impact on MRSI spectral quality.
The current study used a long-term averaging mode in which all rosette rings are traversed before repeating signal averages. In this case, the effects of frequency drift would be more pronounced along the averages dimension rather than spread throughout k-space. In the case of short-term averaging, in which all averages are collected for a single rosette ring before proceeding to the next ring, the opposite would be the case: Frequency drift would be more pronounced throughout k-space than along the averages' dimension. How exactly drift-related artifacts would manifest differently between these two cases is not well understood and should be the subject of future study, either experimentally or via simulation.
In single-voxel MRS, the performance of spectral registration is dependent on voxel size, with small voxels yielding lower SNR, resulting in poorer performance of the spectral registration algorithm. In contrast, the performance of spectral registration in the current MRSI application is not dependent on the spatial resolution of the reconstructed MRSI image. This is because spectral alignment uses the k = 0 FIDs, which contain signal from the whole excited slice/volume. Because MRSI experiments normally excite a large volume, SNR constraints are not expected to be an issue for most applications of the proposed method. This feature may render this method useful in x-nuclear MRS, in which metabolite signals are weaker. Moreover, because time-domain spectral registration uses the early part of the FID for alignment, it is expected to perform well without prior knowledge of the frequencies of the spectral peaks. Chemical shift displacement errors (CSDEs) occur in PRESS-localized MRSI due to the effects of slice-selective localization pulses on metabolite resonances with differing chemical shifts. This can lead to quantification errors in the regions of the volume of interest affected by CSDEs, particularly the edges of the prelocalized volume of interest. Thus, the potential impact of CSDE should still be considered when interpreting MRSI results.

CONCLUSIONS
As rapid spatial-spectral MRSI sequences enter more mainstream use, and given that these sequences may induce large amounts of frequency drift, controlling frequency drift-related artifacts will be an important consideration in the future. We have demonstrated the application of time-domain spectral registration to retrospectively estimate and correct frequency drift errors in rapid, high-resolution rosette MRSI data. The method yielded meaningful improvements in both spectral LW and SNR. ORCID Sneha Senthil https://orcid.org/0000-0002-3601-9497 Jamie Near https://orcid.org/0000-0003-3516-936X TWITTER Jamie Near @JamieNear