Chemical exchange saturation transfer (CEST) imaging is a specific type of magnetic resonance imaging (MRI) technique that is sensitive to exchangeable protons in exogenous or endogenous compounds. After being selectively saturated by radiofrequency pulses at their particular resonance frequencies, the exchangeable labile protons transfer their saturated nuclear magnetization to the bulk water proton pool through chemical exchange. Thus, the degree of bulk water saturation observed becomes sensitive to the labile proton concentration, chemical exchange rate, relaxation times, irradiation parameters, and various other factors (1).
The sensitivity of CEST imaging to exchangeable protons makes it a potential imaging biomarker for various diseases that alter metabolite concentration, such as proteins/peptides (2) and glutamate (3), or alter microenvironment properties, such as pH (4) and temperature (5), as these parameters can directly alter the chemical exchange rates between labile protons and bulk water protons. As a result, CEST MRI has been used to obtain unique imaging contrast in different pathologies, include ischemic stroke and brain tumors (6).
The conventional approach to analyzing CEST MRI data is to examine the asymmetry in the magnetization transfer ratio measured with respect to the resonance frequency of water protons, or asymmetric magnetization transfer ratio (MTRasym) = S( − Δω)/S0 − S(Δω)/S0, where Δω is the frequency shift between the saturation frequency of interest and the frequency of water proton resonance (or the Z-spectrum offset frequency), S(Δω) is the water signal with the application of off-resonance saturation pulses, and S0 is the reference water signal without the application of saturation pulses. This measure, MTRasym, reduces the direct water saturation and the symmetric macromolecular magnetic transfer effects. However, because calculation of MTRasym requires prior knowledge of water proton resonance frequency, which is directly proportional to the static magnetic field B0, traditional measurements of CEST contrast can be highly sensitive to static magnetic field inhomogeneities (Figure 1).
Various approaches have been used to correct B0 inhomogeneities. For example, one method involves acquiring additional Z-spectral points around the theorized water resonance frequency to find the frequency with the lowest signal amplitude, as this should be at or near the water resonance frequency (4, 7). After estimating this local minimum for every voxel, the rest of the Z-spectrum is adjusted accordingly (Figure 2). Although this method is fast and simple, it has several disadvantages including errors introduced because of discretized estimates of ΔB0, restrictions on the detectable range of B0 inhomogeneities, and ΔB0 estimates are easily affected by signal fluctuation because of noise and other contamination. The frequency resolution of the ΔB0 estimates could be improved with fitting and interpolation of Z-spectra, including high-order polynomial fitting (6, 8), smoothing-splines fitting (9), and multipool Lorentzian line fitting (10). Kim et al. (11) proposed water saturation shift referencing (WASSR) method with a similar principle, while requiring an extra WASSR Z-spectra scan with low B1 amplitude to isolate the effect of direct water saturation.
A second method is to use a more traditional approach to B0 field mapping through estimating phase differences from 2 gradient-echo acquisitions with different echo times, or by using a multiecho acquisition scheme (12). A ΔB0 map can be estimated from the difference in the phase measured between 2 different gradient-echo times: ΔB0 = (ϕTE2 − ϕTE1)/γ(TE2 − TE1) (Figure 3). This method overcomes many of the issues associated with the previous technique, but it also has drawbacks including requiring at least 1 voxel where the exact ΔB0 value is known; restrictions on the detectable range of B0 inhomogeneities as determined by the echo times chosen and Nyquist limitations; and, in many cases, multiecho acquisition may not be reasonable (eg, turbo spin echo techniques), so separate acquisitions may be required, resulting in increased acquisition time and possible errors due to dynamic B0 fluctuations (13).
In this study we introduce and demonstrate the benefits of 2 new methods for B0 correction, specifically designed for fast clinical pH-weighted amine proton CEST echo-planar imaging (7, 14), by using k-means clustering combined with Lorentzian fitting (LF) to obtain better-quality images of MTRasym at 3.0 ppm (amine protons) on clinical 3 T magnetic resonance (MR) scanners. In particular, we will compare intrasubject and intersubject consistency, and tumor region signal-to-noise ratio (SNR) of MTRasym at 3.0 ppm between an iterative downsampling estimation (IDE) (15) and a 4D polynomial fitting technique compared with existing approaches.
CEST images for all experiments were acquired using a spin-and-gradient echo (SAGE) echoplanar imaging (EPI) readout (14) for a total of 25 contiguous slices with a slice thickness of 4 mm and no interslice gap. A radiofrequency saturation pulse train, optimized for pH sensitivity by targeting fast-exchanging amine protons (7), consisting of three 100-millisecond Gaussian pulses with peak amplitude B1 = 6 μT were used to sample a total of 29 Z-spectral points around ±3.0 ppm and 0.0 ppm with respect to water resonance frequency (from −3.5 to −2.5 in intervals of 0.1, from −0.3 to +0.3 in intervals of 0.1, and from +2.5 to +3.5 in intervals of 0.1). An additional reference S0 scan with identical parameters and no saturation pulse was acquired with 4 averages. All MRI images were acquired on 3 T scanner (3T PRISMA or SKYRA, Siemens, Erlangen, Germany). The total acquisition time for whole-brain CEST-SAGE-EPI and S0 images was 7 minutes and 30 seconds. In addition to CEST-SAGE-EPI acquisition before contrast administration, all patients received anatomic imaging including T2-weighted fluid-attenuated inversion recovery (FLAIR) images, isotropic 3D T1-weighted scans before and after injection of 0.01 mg/kg Gd-DTPA according to the international standardized brain tumor imaging protocol (16).
A total of 27 adult primary brain tumor patients (33 scans) with high-quality amine CEST-SAGE-EPI and anatomic MRI available as part of standard of care between October 2017 and March 2018 were included in the current study. All patients included in the current study provided IRB-approved informed written consent to have advanced imaging including amine CEST MRI obtained clinically and to allow access to their clinical medical record. Detailed patient characteristics are outlined in Table 1.
ΔB0 Map Estimation
All CEST data were preprocessed with motion correction (mcflirt function, FSL tools) (17).
Central Spectrum Point Correction Method.
Raw CEST data were normalized by the S0 reference scan data, resulting in Z-spectra data. The Z-spectral points acquired were densely sampled around offset frequency of interest (±3.0 ppm for amine protons) and water saturation frequency (0.0 ppm). For every voxel, the signal intensities of the Z-spectral points around the theorized water resonance frequency were compared. The spectral point corresponding to the minimum signal intensity was considered as the B0 offset, and the entire Z-spectrum was adjusted accordingly (Figure 2A).
Z-Spectrum Fitting Correction Method.
Similar to the central spectrum point correction (CSC) method, raw CEST data were first normalized by the S0 reference scan, resulting in Z-spectra data. Then, for every voxel, the Z-spectrum was fitted using a ninth-order polynomial, a smoothing spline, or a single-pool Lorentzian line fit, and the spectrum was interpolated with a frequency sampling step of 0.002 ppm. No parameter constraints were specified for polynomial fitting and smoothing spline fitting (6, 9). The single-pool LF was applied with a constraint of −1 ppm < ΔB0 < +1 ppm. After interpolation, the B0 offset was estimated as the spectral point corresponding to the minimum signal intensity. The Z-spectrum was then adjusted accordingly (Figure 2, C–E).
Multiecho-Phase ΔB0 Mapping.
B0 correction was performed by creating a voxel-wise ΔB0 map from phase data obtained from the 2 gradient echoes in the S0 images, using the following formula: ΔB0 = (phase(TE2) − phase(TE1)) / γ(TE2 − TE1), where TE1 and TE2 are the echo times of the 2 multiecho images and γ is the gyromagnetic ratio (Figure 2B).
IDE Using LF.
A slice-by-slice ΔB0 map estimation was achieved iteratively m times through the following steps (Figure 3A).
Raw CEST data at each offset frequency were normalized with respect to the S0 reference scan data, resulting in Z-spectral measurements.
k-means clustering (kmeans function, MATLAB release 2017b) was used to cluster the pixels with similar patterns in the Z-spectrum. No spatial continuity restriction was applied on the clustering. The number of clusters for each iteration was set to 3m, where m is the iteration number. The original resolution of images can be reached within 7 iterations.
Mean Z-spectra data were fit to the Lorentzian function: c − c / [(x − a)2 + b], where a indicates the shift in the Z-spectrum center point, b relates to the width of the Lorentzian shape, c and b together determine the y-scaling of the function. The fitting started with loosely constrained bounds and narrowed the constraint boundaries after each iteration, m. The constraint bound for ΔB0 estimation was set to be ±1/m ppm. If the goodness of fit was low, R2 < 0.7, then ΔB0 estimation was performed with the conventional CSC method to avoid unrealistic fitting to poor data.
The B0 inhomogeneity (ΔB0) estimation (value of a in the Lorentzian fit) was stored for all pixels, then the Z-spectra data were updated using the newly estimated ΔB0 map.
Determine whether to terminate the iterative process by determining whether the mean square ΔB0 estimates for the current iteration was <10% of the sum of the mean square ΔB0 estimates over all past iterations, indicating little improvement in estimation by the current iteration. (In practice, the iterative ΔB0 estimation typically terminates after 4 or 5 iterations). If this stopping criterion was met, the sum of ΔB0 maps from all iterations was determined to be the final ΔB0 map. If this criterion was not met, steps 2 through 6 were repeated while the number of k-mean clusters were increased and the constraint bounds for ΔB0 estimation were tightened.
Lastly, a 3D Gaussian filter (standard deviation = 1) was applied to the final ΔB0 map to remove potential erroneous estimations.
Two-Stage Lorentzian Estimation With 4D Polynomial Fitting (LF-Poly).
A whole-brain 3D ΔB0 map estimation was achieved through the following 2-stage procedure (Figure 3B).
Raw CEST data were normalized with respect to the S0 reference scan data to calculate the Z-spectra. k-means clustering was used to cluster voxels according to similar patterns in the Z-spectrum. The average Z-spectrum was calculated for all voxels within the same cluster. Lastly, the mean Z-spectrum data were fit with Lorentzian function as described in the previous section.
4D Polynomial Fitting.
ΔB0 estimation from stage 1 was fit to a spatial 4D polynomial function (x, y, z spatial coordinates plus the Z-spectral dimension). This approach was based on common knowledge that ΔB0 is not discontinuous and thus is spatially smooth across the space.
The Z-spectral data were updated using the new estimated ΔB0 map following 4D polynomial fitting. Stage 1 was repeated to correct for the residue of the ΔB0 map. Lastly, the final ΔB0 map estimations was the sum of ΔB0 estimations from the 2 stages, and a 3D Gaussian filter (standard deviation = 1) was then applied to remove potential erroneous estimations.
To optimize the number of k-means clusters and the order of 4D polynomial fitting in the LF-Poly method, combinations of 8 cluster numbers (10, 20, 50, 100, 200, 300, 400, 500) and 8 polynomial orders (5–12) were used to generate ΔB0 estimations and resulting MTRasym images for 7 patient data (the first 7 patients enrolled in the study). The median image quality metrics of the 7 patients' images were then calculated as described in the next section. The performance and efficiency were evaluated to find the optimal combination of parameters.
All methods were compared with the following 2 image quality metrics of ΔB0 map and resulting MTRasym map (at 3.0 ppm): peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM). PSNR was calculated as follows: PSNR = 10log10(peakval2/MSE), where peakval is set to 1 for ΔB0 map and 10 for the MTRasym map, and MSE is the mean square error with the reference image. The ΔB0 maps and resulting MTRasym images generated by the multiecho-phase method were used as ground truth for both image quality metrics. The computation time for each method was also compared. All image processing and image quality metrics' calculation was carried out with MATLAB (Release 2017b, MathWorks, Natick, MA).
The following methods were compared: CSC method; voxel-wise Z-spectrum fitting using ninth-order polynomial, smoothing spline, and single-pool Lorentzian line fit; IDE-LF; and 2-stage LF-Poly.
Medians and standard deviations of the whole-brain MTRasym (at 3.0 ppm) histogram and normal-appearing white matter (NAWM) MTRasym (at 3.0 ppm) of all 27 patients (33 scans) were compared among CSC method, IDE-LF method, multiecho-phase method, and LF-Poly method, by means of Student t test. The MTRasym (at 3.0 ppm) statistics in FLAIR hyperintensity regions were also compared among the aforementioned methods. The SNR of MTRasym in the FLAIR hyperintensity regions was calculated as the ratio between the median MTRasym in FLAIR lesion and the standard deviation of MTRasym in the NAWM region, with Student t test. The NAWM regions of interest (ROIs) and FLAIR hyperintensity ROIs were drawn manually with AFNI software (intensity threshold subsequently applied to segment FLAIR hyperintensity ROIs), and the statistics were calculated with MATLAB (Release 2017b).
An example of the IDE-LF ΔB0 map estimation is shown in Figure 4. This figure illustrates the ΔB0 map estimation from each iteration (Figure 4A) and the sum of ΔB0 map estimation after each iteration (Figure 4B). Note that the number of k-means clusters increased with iteration, and as such, a finer division of brain tissue according to Z-spectra data can be observed. In addition, note that the gradient of ΔB0 also becomes smoother in later iterations. An example of mean square (MS) ΔB0 map – iteration number relationship from one slice is plotted in Figure 4. The MS value of ΔB0 map from each iteration decreases with increase in the iteration number, indicating less extra information about ΔB0 map from later iterations. The stopping criterion of MS(ΔB0) for the current iteration being smaller than 10% of the sum of MS(ΔB0) over all past iterations is met after the fifth iteration in this example. For all 33 scans from 27 patients, the stopping criterion was met after 4.98 (±0.80) iterations (mean ± standard deviation). For the 12 presurgery scans and the 13 postsurgery scans of glioma patients, the stopping criterion was met after 4.91 (±0.73) and 4.82 (±0.81) iterations, respectively. No significant difference was found between the 2 groups (P = .76, Student t test).
The number of k-means clusters and the order of 4D polynomial fitting varied in the range from 10 to 500 and 5 to 12, respectively, to optimize the performance and efficiency of ΔB0 estimation. The PSNR and SSIM metrics of ΔB0 maps and MTRasym images generated with different combinations of parameters are compared in Figure 5. The average of the 4 metrics and the computation time required for each combination of parameters were also plotted in Figure 5. All image quality metrics increased with the increasing number of k-means clusters and the increasing polynomial fitting order. However, the effects of k-means cluster number appeared more influential than those of the polynomial fitting order, and in most circumstances, this increase plateaued when the number of k-means clusters reached between 100 and 200. Although higher values of image quality metrics can be achieved with larger cluster numbers and higher order fitting, the computation time also increased with the increasing number of parameters (Figure 5F). Thus, 2 combinations of parameters were selected for further comparison with the other methods: 100 k-means clusters and eighth-order 4D polynomial fitting and 500 k-means clusters and 12th-order 4D polynomial fitting, as indicated by blue and red circles in Figure 5. The latter corresponds to the combination of parameters with the best image quality performance. The former corresponds to the combination with nearly the same quality (average image quality metrics of 97.97% compared with 99.83%) while requiring less than half of the computation time (26.70 seconds compared with 62.23 seconds per patient).
An example slice of estimated ΔB0 map and resulting MTRasym image generated by different correction methods is illustrated in Figure 6. The quantitative assessments of ΔB0 estimation methods for 7 patients are plotted in Figure 7, including image quality metrics PSNR and SSIM of ΔB0 maps and MTRasym images, computation time, and standard deviation of NAWM MTRasym. The image quality metrics and the standard deviation of NAWM MTRasym show consistent results. Among all methods except for multiecho-phase method, the proposed IDE-LF method and the LF-Poly method perform the best in terms of having the highest PSNR and SSIM values and lowest NAWM MTRasym standard deviations, while the ninth-order polynomial and the smoothing-spline fittings of Z-spectra have the worst performance. Consistent with the quantitative metrics, the proposed IDE-LF method and the LF-Poly method generated MTRasym images that are visually less noisy in Figure 6, while the ninth-order polynomial and the smoothing-spline fitting methods generated highly noisy MTRasym images, particularly around the nasal sinuses, indicating very inaccurate estimation of ΔB0.
The computation time for whole-brain ΔB0 estimation using different methods varies with orders of magnitude. The CSC method and the multiecho-phase method have near real-time execution time. The LF-Poly method has the computation time of ∼≤1 minute. The IDE-LF method, as well as the ninth-order polynomial and the smoothing-spline fitting methods, requires a computation time of 10–15 minutes. Voxel-by-voxel single-pool Lorentzian line fitting requires the longest computation time of ∼40 minutes to 1 hour per patient.
The CSC method, the multiecho-phase method, and the proposed IDE-LF and LF-Poly methods were selected for further comparison of the intertissue (NAWM) variability and intersubject variability, using data of all 27 patients (33 scans). The results are shown in Figure 8, A and B. Among the methods compared, the proposed IDE-LF method and the LF-Poly method show significantly higher intertissue consistency (lower standard deviation of NAWM MTRasym at 3.0 ppm) than the conventional CSC method and the multiecho-phase method (IDE-LF: P = 2.03 × 10−18 compared with the CSC method and P = 1.10 × 10−4 compared with the the multiecho-phase method; LF-Poly(a): P = 9.02 × 10−19 compared with the CSC method and P = 3.03 × 10−4 compared with the multiecho-phase method). There is no significant difference in the standard deviation of NAWM MTRasym at 3.0 ppm between IDE-LF method and the LF-Poly methods (P = .402/.939). The proposed methods also show higher interpatient consistency (lower standard deviation of median NAWM MTRasym at 3.0 ppm) than the CSC and the multiecho-phase methods. The standard deviation of median NAWM MTRasym is 0.2953% for the IDE-LF method and 0.2812%/0.2666% for the LF-Poly methods, compared with 0.3025% for the CSC method, and 2.0011% for the multiecho-phase method. Between the LF-Poly methods with different parameter choices, a slightly lower standard deviation of NAWM MTRasym and a slightly lower standard deviation of median NAWM MTRasym can be observed in the case of 500 k-means clusters and 12th-order polynomial fitting compared with 100 k-means clusters and eighth-order polynomial fitting, but the difference is not significant (P = .2390). A more detailed comparison of NAWM MTRasym statistics can be found in Table 2.
|Method Comparison(NAWM MTRasym SD)||CSCP-value||Multi-echo PhaseP-value||IDE-LFP-value||LF-Poly(a)P-value||LF-Poly(b)P-value|
|Multiecho Phase||P = 6.59 × 10−12***||N/A|
|IDE-LF||P = 2.02 × 10−18***||P = 1.10 × 10−4***||N/A|
|LF-Poly(a)||P = 9.01 × 10−19***||P = 3.03 × 10−4***||P = 0.40||N/A|
|LF-Poly(b)||P = 5.85 × 10−20***||P = 7.49 × 10−5***||P = 0.93||P = 0.23||N/A|
|Median MTRasym at 3 ppm mean (SD)||0.84% (0.30%)||1.34% (2.00%)||0.64% (0.29%)||0.93% (0.28%)||0.88% (0.26%)|
|SD MTRasym at 3 ppm mean (SD)||0.69% (0.12%)||0.51% (0.13%)||0.45% (0.12%)||0.46% (0.12%)||0.45% (0.12%)|
In Figure 9, A and B, the FLAIR hyperintensity lesion MTRasym (at 3.0 ppm) with different B0 correction methods are shown. The IDE-LF method and the LF-Poly methods show significantly higher FLAIR lesion MTRasym SNR than the conventional CSC method (IDE-LF: P = 2.74 × 10−4 compared with the CSC method; LF-kPoly(a): P = 2.08 × 10−5 compared with the CSC method; LF-kPoly(b): P = 3.04 × 10−5 compared with the CSC method). Although the multiecho-phase method also shows higher SNR than the CSC method (not significant, P = .0918), the intersubject deviation of SNR is large, similar to the results in NAWM regions. A more detailed comparison of FLAIR lesion MTRasym statistics can be found in Table 3. With the proposed IDE-LF method, a significant difference in median MTRasym values can be observed in different ROIs and different groups of patients. Presurgery glioma FLAIR lesions (12 scans) have significantly higher MTRasym median than NAWM regions (all scans) with P = 2.35 × 10−10 and postsurgery FLAIR lesions (15 scans) with P = 3.31 × 10−5. Postsurgery FLAIR lesions have significantly lower MTRasym median than presurgery glioma FLAIR lesions (P = .0022).
|Method Comparison(FLAIR Lesion SNR)||CSCP-value||Multiecho PhaseP-value||IDE-LFP-value||LF-Poly(a)P-value||LF-Poly(b)P-value|
|Multiecho Phase||P = 0.0918||N/A|
|IDE-LF||P = 2.7357 × 10−4***||P = 0.2659||N/A|
|LF-Poly(a)||P = 2.0835 × 10−5***||P = 0.2121||P = 0.3365||N/A|
|LF-Poly(b)||P = 3.0385 × 10−5***||P = 0.2389||P = 0.9502||P = 0.0811||N/A|
|FLAIR Lesion SNR mean (SD)||3.1847 (0.8801)||7.2560 (7.1975)||4.4589 (1.5416)||4.2753 (1.1249)||4.4470 (1.3019)|
The computation time of the selected methods (Figure 8C) was consistent with the preliminary results from 7 patient data. The IDE-LF method requires the longest computation time (680.4 seconds ± 342.2 seconds per patient), while the LF-Poly method can provide whole-brain ΔB0 estimation in ∼1 min, 28.6 seconds ± 3.3 seconds for LF-Poly(a) and 65.9 seconds ± 5.3 seconds for LF-Poly(b). The CSC method and the multiecho-phase method have near real-time computation time.
In Figure 8D, the histograms of whole-brain MTRasym using different B0 correction methods were plotted for 4 patients as examples. Compared with the CSC method, the multiecho-phase method, the IDE-LF method, and the LF-Poly method all generated narrower distribution of MTRasym (at 3.0 ppm), possibly indicating improved intertissue consistency. However, the multiecho-phase method shifted the distribution noticeably in some of the patients. The proposed IDE-LF and LF-Poly methods shifted the distribution slightly toward the negative direction, which could be explained by the better ΔB0 estimation close to the sinus regions, alleviating the overestimation of MTRasym (at 3.0 ppm) in areas with very high ΔB0 (Figure 6). No difference of whole-brain MTRasym distribution has been observed between the IDE-LF and LF-Poly methods.
An example of whole-brain ΔB0 correction with multiecho-phase method and IDE-LF method can be seen in Figure 10 to better understand the deviation of ΔB0 estimation between the 2 methods. Several differences have been observed:
The proposed IDE-LF method avoided overestimation of ΔB0 in slice 16 and slice 22, indicated by the more consistent white matter MTRasym across the slices.
The IDE-LF method does not suffer from phase-unwrapping errors, which can be observed in slice 4 for the multiecho-phase method, indicated by the abnormal low MTRasym signal in the right temporal lobe.
This study explored the application of 2 newly proposed methods using k-means clustering combined with LF in B0 correction for clinical CEST images. The first method incorporated iterative downsampling into this estimation, while the second method involved a 4D polynomial fitting to better estimate ΔB0. The results clearly show that the proposed methods outperform the existing methods (CSC method, polynomial/smoothing-spline/single-pool LF, multiecho-phase method) within the context of clinical CEST data that are often undersampled along the Z-spectrum. The proposed methods can serve as promising new approaches for improving CEST image B0 correction, particularly when multiecho CEST images are not available.
Owing to the limit of scan time in clinical settings, clinical CEST data usually suffer from lower SNR and few Z-spectra spectral points compared with research scans, where more number of averages and spectral points can be acquired. Even with accelerated CEST sequences, such as CEST-EPI-SAGE (14), only one measurement of 29 Z-spectral points was possible within clinically allowed time (<10 minute) for a whole-brain-coverage CEST imaging. This limitation restricted the use of many B0 correction methods in clinical use, including WASSR method, which requires excessive scan time.
The low SNR and limited number of spectral points also impair the performance of many local fitting methods, including high-order polynomial fitting and smoothing-spline fitting of Z-spectra. The poor performance shown in Figure 7 was likely because of the error induced in fitting noisy data and unevenly sampled Z-spectral points (Figure 2). The uneven sampling of the Z-spectrum particularly impairs the fitting accuracy of the polynomial method when B0 inhomogeneity is large (outside the range of central spectral points). This can be observed in Figure 6. Although voxel-by-voxel single-pool LF performs better, as it uses the global information of Z-spectra; this method requires unrealistic computation time for whole-brain analysis (40 minutes to 1 hour per patient) and may still suffer from corrupt Z-spectral data.
Compared with the conventional CSC method and the multiecho-phase method, the 2 newly proposed approaches have the advantages of both methods, while avoiding many of their drawbacks. For example, the proposed IDE-LF and LF-Poly methods use the acquired CEST spectra data to perform inline correction, avoiding the potential issue with ΔB0 map fluctuations and added acquisition time. These methods also avoid potential phase-unwrapping errors and inaccuracies in estimation of phase in the MR data. On the other hand, the proposed approaches managed to estimate a wider range of B0 inhomogeneities than the conventional CSC method, as LF uses the shape information of the whole spectrum instead of only central points. The proposed methods also improved the smoothness of the resulting ΔB0 map through an iterative downsampling method in IDE-LF or by 4D spatial polynomial fitting of ΔB0 map in LF-Poly, avoiding potential erroneous values that may arise from other methods.
The improved performance of ΔB0 estimation using k-means clustering combined with LF could be explained by the “denoising” effect on the resulting Z-spectrum data. k-means Clustering of 4D data cluster voxels with similar patterns and features in the Z-spectral data and thus often group similar tissue types together. The averaging of Z-spectral data within the same cluster greatly improves the SNR of the resulting Z-spectral data, thus allowing better estimation of ΔB0 and faster convergence of LF.
In the proposed methods, a single Lorentzian line shape was used to fit the underlying Z-spectral data. This may have caused potential issues with overestimating ΔB0 in voxels with lower values of positive Z-spectrum points (ie, higher MTRasym). The weighting of central spectrum points (from −0.3 ppm to +0.3 ppm) was increased to mitigate this issue, but no overestimation of ΔB0 was observed in any of the patient cases examined. However, we should still keep in mind that the proposed methods may introduce error when applying to CEST data with asymmetric Z-spectra, particularly when the asymmetry approaches free water proton resonance, for example, when investigating the hydroxyl (-OH) proton pool.
The two proposed methods for ΔB0 estimation, IDE-LF and LF-Poly, show similar performance resulting in improved image quality in MTRasym, reduced intra-subject and intersubject NAWM MTRasym deviations, and improved SNR in FLAIR lesion MTRasym. The LF-Poly accounts for the spatial smoothness of ΔB0 inhomogeneities with the 4D polynomial fitting, which at the same time could cause less accurate estimations of local ΔB0 details. However, the much shorter computation time of LF-Poly (∼30 seconds) compared with IDE-LF (∼10 minutes) makes it a better choice for fast generation of MTRasym images in clinical use, while IDE-LF may be considered for a more in-depth quantitative analysis of CEST data when long reconstruction and postprocessing times are not an issue.
No significant difference was observed between the performance of LF-Poly methods with the 2 different parameter settings (100 k-means clusters + eighth-order polynomial fitting and 500 k-means clusters + 12th-order polynomial fitting). The plateau of performance with increasing number of clusters and polynomial fitting orders could be explained by the failure to converge within limited iterations for large number of clusters and overfitting of the model with singular covariance matrices for limited number of Z-spectral points. The dependency of performance on the number of clusters and polynomial fitting orders is similar in all 7 patients' data used for optimization (data not shown), despite the different tumor types, sizes, shapes, locations, grades, and treatment status. This indicates that the choice of optimal parameters would be minimally affected by the presence and characteristics of the brain tumors. The parameters of 100 k-means clusters and eighth-order polynomial fitting require less computation time and thus fit better the clinical need.
One limitation of this study is that we validated only our method with one specific CEST MRI protocol and limited patient cohort (brain tumor patients). Different CEST parameters, including Z-spectral sampling, acquisition matrix size, and number of averages, could potentially change the performance and optimization of parameters. The optimized parameters would also be expected to be different in phantom experiments and in different organs. However, once the optimized parameters are determined for a particular CEST protocol, the parameters should be applicable to any subject. Another limitation is the small number of patient samples examined. With the small sample size, we were not able to give a conclusion on the dependency of performance, parameter optimization, and stopping criteria, on the size and location of brain tumors.
The proposed methods could be further improved by optimizing the parameter choices (number of clusters for each iteration, stopping criteria, constraint bound for LF, etc.) to achieve higher efficiency for high-quality ΔB0 map generation. It would also be worthwhile to explore the possibility of reducing the number of central spectrum points needed, while maintaining the reliability of generated ΔB0 map.
In conclusion, we proposed 2 new methods for CEST B0 correction using k-means clustering combined with LF for improving B0 correction for fast clinical pH-weighted CEST imaging: the iterative downsampling estimation using Lorentzian fitting (IDE-LF) and the 2-stage Lorentzian estimation with 4D polynomial fitting (LF-Poly). The proposed methods were shown to improve the image quality, intertissue consistency, interpatient consistency, and lesion SNR of MTRasym (at 3.0 ppm) images.