Introduction
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 (MTR_{asym}) = S( − Δω)/S_{0} − S(Δω)/S_{0}, where Δω is the frequency shift between the saturation frequency of interest and the frequency of water proton resonance (or the Zspectrum offset frequency), S(Δω) is the water signal with the application of offresonance saturation pulses, and S_{0} is the reference water signal without the application of saturation pulses. This measure, MTR_{asym}, reduces the direct water saturation and the symmetric macromolecular magnetic transfer effects. However, because calculation of MTR_{asym} requires prior knowledge of water proton resonance frequency, which is directly proportional to the static magnetic field B_{0}, traditional measurements of CEST contrast can be highly sensitive to static magnetic field inhomogeneities (Figure 1).
Figure 1.
Example of asymmetric magnetization transfer ratio (MTR_{asym} at 3.0 ppm) map with and without B_{0} correction. (A) and (C) demonstrate MTR_{asym} (at 3.0 ppm) of the same image slice with and without B_{0} correction, respectively. (B) is the ΔB_{0} map used for correction. Note the higher extent of MTR_{asym} contrast in (A) compared with that after B_{0} correction in (C).
Various approaches have been used to correct B_{0} inhomogeneities. For example, one method involves acquiring additional Zspectral 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 Zspectrum is adjusted accordingly (Figure 2). Although this method is fast and simple, it has several disadvantages including errors introduced because of discretized estimates of ΔB_{0}, restrictions on the detectable range of B_{0} inhomogeneities, and ΔB_{0} estimates are easily affected by signal fluctuation because of noise and other contamination. The frequency resolution of the ΔB_{0} estimates could be improved with fitting and interpolation of Zspectra, including highorder polynomial fitting (6, 8), smoothingsplines 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 Zspectra scan with low B1 amplitude to isolate the effect of direct water saturation.
Figure 2.
Example of Zspectrum B_{0} correction with different methods. Central spectral points correction (CSC) methods (A): The minimum point of the original Zspectrum (blue) is found to be at +0.2 ppm. The whole Zspectrum is shifted so that the minimum point occurs at 0 ppm. The ΔB_{0} map generated by this method is shown on the right. Multiechophase method (B): the ΔB_{0} map is generated from the subtraction of phase images at 2 different echo times, TE1 and TE2. Highorder polynomial (C), smoothingspline (D), singlepool Lorenztian fitting (LF) method (E): the Zspectrum data points are fitted with a ninthorder polynomial function (C)/smoothingspline function (D)/singlepool Lorenztian function (E). The minimum points of the fitting curves are considered to be the B_{0} offsets.
A second method is to use a more traditional approach to B_{0} field mapping through estimating phase differences from 2 gradientecho acquisitions with different echo times, or by using a multiecho acquisition scheme (12). A ΔB_{0} map can be estimated from the difference in the phase measured between 2 different gradientecho times: ΔB_{0} = (ϕ_{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 ΔB_{0} value is known; restrictions on the detectable range of B_{0} 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 B_{0} fluctuations (13).
Figure 3.
Flowchart of the proposed ΔB_{0} estimation methods. The CEST data were preprocessed with motion correction and normalization to reference scans. The normalized CEST data were then used to perform ΔB_{0} estimation with the 2 newly proposed methods. The details of the methods are described in the Methodology.
In this study we introduce and demonstrate the benefits of 2 new methods for B_{0} correction, specifically designed for fast clinical pHweighted amine proton CEST echoplanar imaging (7, 14), by using kmeans clustering combined with Lorentzian fitting (LF) to obtain betterquality images of MTR_{asym} 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 signaltonoise ratio (SNR) of MTR_{asym} at 3.0 ppm between an iterative downsampling estimation (IDE) (15) and a 4D polynomial fitting technique compared with existing approaches.
Methodology
MRI Acquisition
CEST images for all experiments were acquired using a spinandgradient 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 fastexchanging amine protons (7), consisting of three 100millisecond Gaussian pulses with peak amplitude B_{1} = 6 μT were used to sample a total of 29 Zspectral 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 S_{0} 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 wholebrain CESTSAGEEPI and S_{0} images was 7 minutes and 30 seconds. In addition to CESTSAGEEPI acquisition before contrast administration, all patients received anatomic imaging including T2weighted fluidattenuated inversion recovery (FLAIR) images, isotropic 3D T1weighted scans before and after injection of 0.01 mg/kg GdDTPA according to the international standardized brain tumor imaging protocol (16).
Patients
A total of 27 adult primary brain tumor patients (33 scans) with highquality amine CESTSAGEEPI 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 IRBapproved 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.
Table 1.
Characteristics 
All Patients (N = 27) 
Age, Year 

Median (Interquartile Range) 
52 (47–61.5) 
Range 
19–76 
Sex, No. (%) 

Male 
13 (48.1%) 
Female 
14 (51.9%) 
Diagnosis, No. (%) 

Glioblastoma (WHO Grade IV Glioma) 
10 (37.0%) 
Meningioma 
7 (25.9%) 
WHO Grade II Glioma 
6 (22.2%) 
WHO Grade III Glioma 
3 (11.1%) 
Gliosarcoma 
1 (3.7%) 
ΔB_{0} 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 S_{0} reference scan data, resulting in Zspectra data. The Zspectral 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 Zspectral points around the theorized water resonance frequency were compared. The spectral point corresponding to the minimum signal intensity was considered as the B_{0} offset, and the entire Zspectrum was adjusted accordingly (Figure 2A).
ZSpectrum Fitting Correction Method.
Similar to the central spectrum point correction (CSC) method, raw CEST data were first normalized by the S_{0} reference scan, resulting in Zspectra data. Then, for every voxel, the Zspectrum was fitted using a ninthorder polynomial, a smoothing spline, or a singlepool 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 singlepool LF was applied with a constraint of −1 ppm < ΔB_{0} < +1 ppm. After interpolation, the B_{0} offset was estimated as the spectral point corresponding to the minimum signal intensity. The Zspectrum was then adjusted accordingly (Figure 2, C–E).
MultiechoPhase ΔB_{0} Mapping.
B_{0} correction was performed by creating a voxelwise ΔB_{0} map from phase data obtained from the 2 gradient echoes in the S_{0} images, using the following formula: ΔB_{0} = (phase(TE_{2}) − phase(TE_{1})) / γ(TE_{2} − TE_{1}), where TE_{1} and TE_{2} are the echo times of the 2 multiecho images and γ is the gyromagnetic ratio (Figure 2B).
IDE Using LF.
A slicebyslice ΔB_{0} map estimation was achieved iteratively m times through the following steps (Figure 3A).
Step 1.
Raw CEST data at each offset frequency were normalized with respect to the S_{0} reference scan data, resulting in Zspectral measurements.
Step 2.
kmeans clustering (kmeans function, MATLAB release 2017b) was used to cluster the pixels with similar patterns in the Zspectrum. No spatial continuity restriction was applied on the clustering. The number of clusters for each iteration was set to 3^{m}, where m is the iteration number. The original resolution of images can be reached within 7 iterations.
Step 3.
The average Zspectra data were averaged across all pixels within the same cluster.
Step 4.
Mean Zspectra data were fit to the Lorentzian function: c − c / [(x − a)^{2} + b], where a indicates the shift in the Zspectrum center point, b relates to the width of the Lorentzian shape, c and b together determine the yscaling of the function. The fitting started with loosely constrained bounds and narrowed the constraint boundaries after each iteration, m. The constraint bound for ΔB_{0} estimation was set to be ±1/m ppm. If the goodness of fit was low, R^{2} < 0.7, then ΔB_{0} estimation was performed with the conventional CSC method to avoid unrealistic fitting to poor data.
Step 5.
The B_{0} inhomogeneity (ΔB_{0}) estimation (value of a in the Lorentzian fit) was stored for all pixels, then the Zspectra data were updated using the newly estimated ΔB_{0} map.
Step 6.
Determine whether to terminate the iterative process by determining whether the mean square ΔB_{0} estimates for the current iteration was <10% of the sum of the mean square ΔB_{0} estimates over all past iterations, indicating little improvement in estimation by the current iteration. (In practice, the iterative ΔB_{0} estimation typically terminates after 4 or 5 iterations). If this stopping criterion was met, the sum of ΔB_{0} maps from all iterations was determined to be the final ΔB_{0} map. If this criterion was not met, steps 2 through 6 were repeated while the number of kmean clusters were increased and the constraint bounds for ΔB_{0} estimation were tightened.
Step 7.
Lastly, a 3D Gaussian filter (standard deviation = 1) was applied to the final ΔB_{0} map to remove potential erroneous estimations.
TwoStage Lorentzian Estimation With 4D Polynomial Fitting (LFPoly).
A wholebrain 3D ΔB_{0} map estimation was achieved through the following 2stage procedure (Figure 3B).
Stage 1.
Raw CEST data were normalized with respect to the S_{0} reference scan data to calculate the Zspectra. kmeans clustering was used to cluster voxels according to similar patterns in the Zspectrum. The average Zspectrum was calculated for all voxels within the same cluster. Lastly, the mean Zspectrum data were fit with Lorentzian function as described in the previous section.
4D Polynomial Fitting.
ΔB_{0} estimation from stage 1 was fit to a spatial 4D polynomial function (x, y, z spatial coordinates plus the Zspectral dimension). This approach was based on common knowledge that ΔB_{0} is not discontinuous and thus is spatially smooth across the space.
Stage 2.
The Zspectral data were updated using the new estimated ΔB_{0} map following 4D polynomial fitting. Stage 1 was repeated to correct for the residue of the ΔB_{0} map. Lastly, the final ΔB_{0} map estimations was the sum of ΔB_{0} 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 kmeans clusters and the order of 4D polynomial fitting in the LFPoly method, combinations of 8 cluster numbers (10, 20, 50, 100, 200, 300, 400, 500) and 8 polynomial orders (5–12) were used to generate ΔB_{0} estimations and resulting MTR_{asym} 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.
Methods' Comparison
All methods were compared with the following 2 image quality metrics of ΔB_{0} map and resulting MTR_{asym} map (at 3.0 ppm): peak signaltonoise ratio (PSNR) and structural similarity index (SSIM). PSNR was calculated as follows: PSNR = 10log_{10}(peakval^{2}/MSE), where peakval is set to 1 for ΔB_{0} map and 10 for the MTR_{asym} map, and MSE is the mean square error with the reference image. The ΔB_{0} maps and resulting MTR_{asym} images generated by the multiechophase 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; voxelwise Zspectrum fitting using ninthorder polynomial, smoothing spline, and singlepool Lorentzian line fit; IDELF; and 2stage LFPoly.
Medians and standard deviations of the wholebrain MTR_{asym} (at 3.0 ppm) histogram and normalappearing white matter (NAWM) MTR_{asym} (at 3.0 ppm) of all 27 patients (33 scans) were compared among CSC method, IDELF method, multiechophase method, and LFPoly method, by means of Student t test. The MTR_{asym} (at 3.0 ppm) statistics in FLAIR hyperintensity regions were also compared among the aforementioned methods. The SNR of MTR_{asym} in the FLAIR hyperintensity regions was calculated as the ratio between the median MTR_{asym} in FLAIR lesion and the standard deviation of MTR_{asym} 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).
Results
An example of the IDELF ΔB_{0} map estimation is shown in Figure 4. This figure illustrates the ΔB_{0} map estimation from each iteration (Figure 4A) and the sum of ΔB_{0} map estimation after each iteration (Figure 4B). Note that the number of kmeans clusters increased with iteration, and as such, a finer division of brain tissue according to Zspectra data can be observed. In addition, note that the gradient of ΔB_{0} also becomes smoother in later iterations. An example of mean square (MS) ΔB_{0} map – iteration number relationship from one slice is plotted in Figure 4. The MS value of ΔB_{0} map from each iteration decreases with increase in the iteration number, indicating less extra information about ΔB_{0} map from later iterations. The stopping criterion of MS(ΔB_{0}) for the current iteration being smaller than 10% of the sum of MS(ΔB_{0}) 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).
Figure 4.
Example of Iterative downsampling estimation (IDELF) of ΔB_{0} map. (A) shows the ΔB_{0} map estimation from iterations 1–7 with kmeans clustering of increasing cluster numbers. The mean square (MS) of the ΔB_{0} map estimation decreases with the increase in iteration number, indicating less extra information of the ΔB_{0} map obtained from each iteration. An example of the relationship between MS values and iteration numbers for one slice of the brain is plotted. (B) demonstrates the sum of ΔB_{0} estimations after iterations 1–7. The ΔB_{0} estimations become smoother and contain more details as the iteration number increases. The last image shows the final ΔB_{0} map after Gaussian filtering.
The number of kmeans 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 ΔB_{0} estimation. The PSNR and SSIM metrics of ΔB_{0} maps and MTR_{asym} 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 kmeans clusters and the increasing polynomial fitting order. However, the effects of kmeans cluster number appeared more influential than those of the polynomial fitting order, and in most circumstances, this increase plateaued when the number of kmeans 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 kmeans clusters and eighthorder 4D polynomial fitting and 500 kmeans clusters and 12thorder 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).
Figure 5.
Image quality metrics and computation time for the 2stage Lorentzian estimation with 4D polynomial fitting (LFPoly) method with differing numbers of kmeans clusters (range from 10 to 500) and orders of 4D polynomial fitting (range from 5 to 12). The figure shows the PSNR (A) and SSIM (B) metrics of ΔB_{0} maps and MTR_{asym} images (C, D) generated with different combinations of parameters. All image quality metrics are normalized to the highest values of all combinations of parameters. The average of the 4 metrics is plotted in (E). The computation time required for each combination of parameters is plotted in (F). The blue and red circles correspond to the cases with 100 kmeans clusters + eighthorder 4D polynomial fitting and 500 kmeans clusters + 12thorder 4D polynomial fitting, respectively. The 2 parameter combinations were used for further comparison with other methods.
An example slice of estimated ΔB_{0} map and resulting MTR_{asym} image generated by different correction methods is illustrated in Figure 6. The quantitative assessments of ΔB_{0} estimation methods for 7 patients are plotted in Figure 7, including image quality metrics PSNR and SSIM of ΔB_{0} maps and MTR_{asym} images, computation time, and standard deviation of NAWM MTR_{asym}. The image quality metrics and the standard deviation of NAWM MTR_{asym} show consistent results. Among all methods except for multiechophase method, the proposed IDELF method and the LFPoly method perform the best in terms of having the highest PSNR and SSIM values and lowest NAWM MTR_{asym} standard deviations, while the ninthorder polynomial and the smoothingspline fittings of Zspectra have the worst performance. Consistent with the quantitative metrics, the proposed IDELF method and the LFPoly method generated MTR_{asym} images that are visually less noisy in Figure 6, while the ninthorder polynomial and the smoothingspline fitting methods generated highly noisy MTR_{asym} images, particularly around the nasal sinuses, indicating very inaccurate estimation of ΔB_{0}.
Figure 6.
ΔB_{0} maps (A) and resulting MTR_{asym} images (B) generated by different methods. IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting. Black arrow: MTR_{asym} (at 3.0 ppm) contrast in tumor area; white arrow: region with severe B_{0} inhomogeneity.
Figure 7.
Image quality metrics for ΔB_{0} maps and resulting MTR_{asym} (at 3.0 ppm) images generated by different methods (data from 7 patients). The PSNR and SSIM metrics of ΔB_{0} maps and MTR_{asym} (at 3.0 ppm) images are shown in (A–D). The computation times of different methods are plotted in (E). (F) shows the standard deviation of normal appearing white matter (NAWM) MTR_{asym} (at 3.0 ppm). Subplots were included in (E) and (F) to zoom in the metric scale for selected methods. Triangle markers of different colors represent data points from different patients. CSC: central spectrum points correction method; IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 kmeans clusters and eighthorder polynomial and (b) parameters of 500 kmeans clusters and 12thorder polynomial. Box plot: the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles.
The computation time for wholebrain ΔB_{0} estimation using different methods varies with orders of magnitude. The CSC method and the multiechophase method have near realtime execution time. The LFPoly method has the computation time of ∼≤1 minute. The IDELF method, as well as the ninthorder polynomial and the smoothingspline fitting methods, requires a computation time of 10–15 minutes. Voxelbyvoxel singlepool Lorentzian line fitting requires the longest computation time of ∼40 minutes to 1 hour per patient.
The CSC method, the multiechophase method, and the proposed IDELF and LFPoly 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 IDELF method and the LFPoly method show significantly higher intertissue consistency (lower standard deviation of NAWM MTR_{asym} at 3.0 ppm) than the conventional CSC method and the multiechophase method (IDELF: P = 2.03 × 10^{−18} compared with the CSC method and P = 1.10 × 10^{−4} compared with the the multiechophase method; LFPoly(a): P = 9.02 × 10^{−19} compared with the CSC method and P = 3.03 × 10^{−4} compared with the multiechophase method). There is no significant difference in the standard deviation of NAWM MTR_{asym} at 3.0 ppm between IDELF method and the LFPoly methods (P = .402/.939). The proposed methods also show higher interpatient consistency (lower standard deviation of median NAWM MTR_{asym} at 3.0 ppm) than the CSC and the multiechophase methods. The standard deviation of median NAWM MTR_{asym} is 0.2953% for the IDELF method and 0.2812%/0.2666% for the LFPoly methods, compared with 0.3025% for the CSC method, and 2.0011% for the multiechophase method. Between the LFPoly methods with different parameter choices, a slightly lower standard deviation of NAWM MTR_{asym} and a slightly lower standard deviation of median NAWM MTR_{asym} can be observed in the case of 500 kmeans clusters and 12thorder polynomial fitting compared with 100 kmeans clusters and eighthorder polynomial fitting, but the difference is not significant (P = .2390). A more detailed comparison of NAWM MTR_{asym} statistics can be found in Table 2.
Figure 8.
MTR_{asym} statistics with selected B_{0} correction methods (data from 27 patients). Standard deviations (A) and medians of NAWM MTR_{asym} (at 3.0 ppm) (B) are plotted. The IDELF method shows the lowest standard deviation of NAWM MTR_{asym} (at 3.0 ppm) and the lowest intersubject deviation of median NAWM MTR_{asym} (at 3.0 ppm) (***P < .001). The computation time of wholebrain ΔB_{0} estimation with different methods is shown in (C). The CSC method and the multiechophase method require the least computation time (median computation time, 1.0 and 0.7 seconds, respectively), while the IDELF method requires the most (median computation time, 697.2 seconds). Examples of wholebrain MTR_{asym} (at 3.0 ppm) histograms of 4 patients' data is demonstrated in (D). Generally speaking, the IDELF method and the LFPoly method lead to narrower distribution, while multiechophase method shifts the distribution in some of the cases. CSC: central spectrum points correction method; IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 kmeans clusters and 8th order polynomial, (b) with parameters of 500 kmeans clusters and 12th order polynomial. Box plot: the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles.
Figure 9.
Comparison of ROI dependence on MTRasym statistics. FLAIR hyperintensity region of interest (ROI) MTR_{asym} (at 3.0 ppm) SNR (calculated as median of MTR_{asym} / standard deviation of NAWM MTR_{asym}) in presurgery scans of patients with glioma is shown in (A). The IDELF method and the LFPoly method shows significantly higher SNR than the conventional CSC method (***P < .001). The phase method also shows higher SNR. However, the intersubject deviation is large, similar to the case in NAWM MTR_{asym}. Comparison of MTR_{asym} (at 3.0 ppm) median values in NAWM, FLAIR lesion in presurgery scans of patients with glioma, FLAIR lesion in postsurgery scans (postOP) of patients with glioma, with IDELF B_{0} correction method is shown in (B). NAWM regions have significantly lower MTR_{asym} median compared with presurgery glioma FLAIR lesions (***P < .001) and with postsurgery FLAIR lesions (***P < .001). Postsurgery FLAIR lesions have significantly lower MTR_{asym} median than presurgery glioma FLAIR lesions (**P < .01). CSC: central spectrum points correction method; IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 kmeans clusters and 8^{th} order polynomial, (b) with parameters of 500 kmeans clusters and 12^{th} order polynomial. Box plot: the central mark indicates the median, and the bottom and top edges of the box indicate the 25^{th} and 75^{th} percentiles.
Table 2.
NAWM MTR_{asym} (at 3.0 ppm) Statistics
Method Comparison(NAWM MTR_{asym} SD) 
CSCPvalue 
Multiecho PhasePvalue 
IDELFPvalue 
LFPoly(a)Pvalue 
LFPoly(b)Pvalue 
CSC 
N/A 




Multiecho Phase 
P = 6.59 × 10^{−12}***

N/A 



IDELF 
P = 2.02 × 10^{−18}***

P = 1.10 × 10^{−4}***

N/A 


LFPoly(a) 
P = 9.01 × 10^{−19}***

P = 3.03 × 10^{−4}***

P = 0.40 
N/A 

LFPoly(b) 
P = 5.85 × 10^{−20}***

P = 7.49 × 10^{−5}***

P = 0.93 
P = 0.23 
N/A 
Median MTR_{asym} 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 MTR_{asym} 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 MTR_{asym} (at 3.0 ppm) with different B_{0} correction methods are shown. The IDELF method and the LFPoly methods show significantly higher FLAIR lesion MTR_{asym} SNR than the conventional CSC method (IDELF: P = 2.74 × 10^{−4} compared with the CSC method; LFkPoly(a): P = 2.08 × 10^{−5} compared with the CSC method; LFkPoly(b): P = 3.04 × 10^{−5} compared with the CSC method). Although the multiechophase 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 MTR_{asym} statistics can be found in Table 3. With the proposed IDELF method, a significant difference in median MTR_{asym} values can be observed in different ROIs and different groups of patients. Presurgery glioma FLAIR lesions (12 scans) have significantly higher MTR_{asym} 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 MTR_{asym} median than presurgery glioma FLAIR lesions (P = .0022).
Table 3.
Presurgery FLAIR Lesion MTR_{asym} (at 3.0 ppm) Statistics in Patients With Glioma
Method Comparison(FLAIR Lesion SNR) 
CSCPvalue 
Multiecho PhasePvalue 
IDELFPvalue 
LFPoly(a)Pvalue 
LFPoly(b)Pvalue 
CSC 
N/A 




Multiecho Phase 
P = 0.0918 
N/A 



IDELF 
P = 2.7357 × 10^{−4}***

P = 0.2659 
N/A 


LFPoly(a) 
P = 2.0835 × 10^{−5}***

P = 0.2121 
P = 0.3365 
N/A 

LFPoly(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 IDELF method requires the longest computation time (680.4 seconds ± 342.2 seconds per patient), while the LFPoly method can provide wholebrain ΔB_{0} estimation in ∼1 min, 28.6 seconds ± 3.3 seconds for LFPoly(a) and 65.9 seconds ± 5.3 seconds for LFPoly(b). The CSC method and the multiechophase method have near realtime computation time.
In Figure 8D, the histograms of wholebrain MTR_{asym} using different B_{0} correction methods were plotted for 4 patients as examples. Compared with the CSC method, the multiechophase method, the IDELF method, and the LFPoly method all generated narrower distribution of MTR_{asym} (at 3.0 ppm), possibly indicating improved intertissue consistency. However, the multiechophase method shifted the distribution noticeably in some of the patients. The proposed IDELF and LFPoly methods shifted the distribution slightly toward the negative direction, which could be explained by the better ΔB_{0} estimation close to the sinus regions, alleviating the overestimation of MTR_{asym} (at 3.0 ppm) in areas with very high ΔB_{0} (Figure 6). No difference of wholebrain MTR_{asym} distribution has been observed between the IDELF and LFPoly methods.
An example of wholebrain ΔB_{0} correction with multiechophase method and IDELF method can be seen in Figure 10 to better understand the deviation of ΔB_{0} estimation between the 2 methods. Several differences have been observed:
The proposed IDELF method avoided overestimation of ΔB_{0} in slice 16 and slice 22, indicated by the more consistent white matter MTR_{asym} across the slices.
The IDELF method does not suffer from phaseunwrapping errors, which can be observed in slice 4 for the multiechophase method, indicated by the abnormal low MTR_{asym} signal in the right temporal lobe.
Figure 10.
Wholebrain ΔB_{0} maps and resulting MTR_{asym} (at 3.0 ppm) maps generated by the multiechophase method (A, B) and the IDELF method (C, D). The MTR_{asym} (at 3.0 ppm) maps generated by the IDELF method show more consistency across the slices compared with the multiechophase method.
Discussion
This study explored the application of 2 newly proposed methods using kmeans clustering combined with LF in B_{0} 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 ΔB_{0}. The results clearly show that the proposed methods outperform the existing methods (CSC method, polynomial/smoothingspline/singlepool LF, multiechophase method) within the context of clinical CEST data that are often undersampled along the Zspectrum. The proposed methods can serve as promising new approaches for improving CEST image B_{0} 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 Zspectra spectral points compared with research scans, where more number of averages and spectral points can be acquired. Even with accelerated CEST sequences, such as CESTEPISAGE (14), only one measurement of 29 Zspectral points was possible within clinically allowed time (<10 minute) for a wholebraincoverage CEST imaging. This limitation restricted the use of many B_{0} 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 highorder polynomial fitting and smoothingspline fitting of Zspectra. The poor performance shown in Figure 7 was likely because of the error induced in fitting noisy data and unevenly sampled Zspectral points (Figure 2). The uneven sampling of the Zspectrum particularly impairs the fitting accuracy of the polynomial method when B_{0} inhomogeneity is large (outside the range of central spectral points). This can be observed in Figure 6. Although voxelbyvoxel singlepool LF performs better, as it uses the global information of Zspectra; this method requires unrealistic computation time for wholebrain analysis (40 minutes to 1 hour per patient) and may still suffer from corrupt Zspectral data.
Compared with the conventional CSC method and the multiechophase method, the 2 newly proposed approaches have the advantages of both methods, while avoiding many of their drawbacks. For example, the proposed IDELF and LFPoly methods use the acquired CEST spectra data to perform inline correction, avoiding the potential issue with ΔB_{0} map fluctuations and added acquisition time. These methods also avoid potential phaseunwrapping 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 B_{0} 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 ΔB_{0} map through an iterative downsampling method in IDELF or by 4D spatial polynomial fitting of ΔB_{0} map in LFPoly, avoiding potential erroneous values that may arise from other methods.
The improved performance of ΔB_{0} estimation using kmeans clustering combined with LF could be explained by the “denoising” effect on the resulting Zspectrum data. kmeans Clustering of 4D data cluster voxels with similar patterns and features in the Zspectral data and thus often group similar tissue types together. The averaging of Zspectral data within the same cluster greatly improves the SNR of the resulting Zspectral data, thus allowing better estimation of ΔB_{0} and faster convergence of LF.
In the proposed methods, a single Lorentzian line shape was used to fit the underlying Zspectral data. This may have caused potential issues with overestimating ΔB_{0} in voxels with lower values of positive Zspectrum points (ie, higher MTR_{asym}). The weighting of central spectrum points (from −0.3 ppm to +0.3 ppm) was increased to mitigate this issue, but no overestimation of ΔB_{0} 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 Zspectra, particularly when the asymmetry approaches free water proton resonance, for example, when investigating the hydroxyl (OH) proton pool.
The two proposed methods for ΔB_{0} estimation, IDELF and LFPoly, show similar performance resulting in improved image quality in MTR_{asym}, reduced intrasubject and intersubject NAWM MTR_{asym} deviations, and improved SNR in FLAIR lesion MTR_{asym}. The LFPoly accounts for the spatial smoothness of ΔB_{0} inhomogeneities with the 4D polynomial fitting, which at the same time could cause less accurate estimations of local ΔB_{0} details. However, the much shorter computation time of LFPoly (∼30 seconds) compared with IDELF (∼10 minutes) makes it a better choice for fast generation of MTR_{asym} images in clinical use, while IDELF may be considered for a more indepth quantitative analysis of CEST data when long reconstruction and postprocessing times are not an issue.
No significant difference was observed between the performance of LFPoly methods with the 2 different parameter settings (100 kmeans clusters + eighthorder polynomial fitting and 500 kmeans clusters + 12^{th}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 Zspectral 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 kmeans clusters and eighthorder 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 Zspectral 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 highquality ΔB_{0} 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 ΔB_{0} map.
In conclusion, we proposed 2 new methods for CEST B_{0} correction using kmeans clustering combined with LF for improving B_{0} correction for fast clinical pHweighted CEST imaging: the iterative downsampling estimation using Lorentzian fitting (IDELF) and the 2stage Lorentzian estimation with 4D polynomial fitting (LFPoly). The proposed methods were shown to improve the image quality, intertissue consistency, interpatient consistency, and lesion SNR of MTR_{asym} (at 3.0 ppm) images.
Acknowledgements
We gratefully acknowledge the support of the American Cancer Society (ACS), Art of the Brain, University of California Research Coordinating Committee, UCLA Jonsson Comprehensive Cancer Center Seed Grant, UCLA SPORE, and NIH/NCI Funding. The study was funded by the American Cancer Society (ACS) Research Scholar Grant (RSG1500301CCE) (Ellingson); Art of the Brain (Cloughesy); University of California Research Coordinating Committee (Ellingson); UCLA Jonsson Comprehensive Cancer Center Seed Grant (Ellingson); UCLA SPORE in Brain Cancer (NIH/NCI 1P50CA21101501A1) (Ellingson, Liau, Nghiemphu, Lai, Pope, Cloughesy); and NIH/NCI 1R21CA22375701 (Ellingson).
Disclosures: No disclosures to report.
Conflict of Interest: The authors have no conflict of interest to declare.
References

Zaiss M, Bachert P. Exchangedependent relaxation in the rotating frame for slow and intermediate exchangemodeling offresonant spinlock and chemical exchange saturation transfer. NMR Biomed. 2013;26:507–518.

Zhou JY, Blakeley JO, Hua J, Kim M, Laterra J, Pomper MG, van Zijl PC. Practical data acquisition method for human brain tumor amide proton transfer (APT) imaging. Magn Reson Med. 2008;60:842–849.

Cai KJ, Haris M, Singh A, Kogan F, Greenberg JH, Hariharan H, Detre JA, Reddy R. Magnetic resonance imaging of glutamate. Nat Med. 2012;18:302–306.

Harris RJ, Cloughesy TF, Liau LM, Prins RM, Antonios JP, Li D, Yong WH, Pope WB, Lai A, Nghiemphu PL, Ellingson BM. pHWeighted molecular imaging of gliomas using amine chemical exchange saturation transfer MRI. Neuro Oncol. 2015;17:1514–1524.

Langereis S, Keupp J, van Velthoven JL, de Roos IH, Burdinski D, Pikkemaat JA, Grüll H. A temperaturesensitive liposomal 1H CEST and 19F contrast agent for MR imageguided drug delivery. J Am Chem Soc. 2009;131:1380–1381.

Zhou J, Lal B, Wilson DA, Laterra J, van Zijl PC. Amide proton transfer (APT) contrast for imaging of brain tumors. Magn Reson Med. 2003;50:1120–1126.

Harris RJ, Cloughesy TF, Liau LM, Nghiemphu PL, Lai A, Pope WB, Ellingson BM. Simulation, phantom validation, and clinical evaluation of fast pHweighted molecular imaging using amine chemical exchange saturation transfer echo planar imaging (CESTEPI) in glioma at 3 T. NMR Biomed. 2016;29:1563–1576.

Zhou J, Payen JF, Wilson DA, Traystman RJ, van Zijl PC. Using the amide proton signals of intracellular proteins and peptides to detect pH effects in MRI. Nat Med. 2003;9:1085–1090.

Stancanello J, Terreno E, Castelli DD, Cabella C, Uggeri F, Aime S. Development and validation of a smoothingsplinesbased correction method for improving the analysis of CESTMR images. Contrast Media Mol Imaging. 2008;3:136–149.

Zaiss M, Schmitt B, Bachert P. Quantitative separation of CEST effect from magnetization transfer and spillover effects by Lorentzianlinefit analysis of zspectra. J Magn Reson. 2011;211:149–155.

Kim M, Gillen J, Landman BA, Zhou J, van Zijl PC. Water saturation shift referencing (WASSR) for chemical exchange saturation transfer (CEST) experiments. Magn Reson Med. 2009;61:1441–1450.

Sun PZ, Farrar CT, Sorensen AG. Correction for artifacts induced by B0 and B1 field inhomogeneities in pHsensitive chemical exchange saturation transfer (CEST) Imaging. Magn Reson Med. 2007;58:1207–1215.

Durand E, van de Moortele PF, PachotClouard M, Le Bihan D. Artifact due to Bo fluctuations in fMRI: correction using the kspace central line. Magn Reson Med. 2001;46:198–201.

Harris RJ, Yao J, Chakhoyan A, Raymond C, Leu K, Liau LM, Nghiemphu PL, Lai A, Salamon N, Pope WB, Cloughesy TF, Ellingson BM. Simultaneous pHsensitive and oxygensensitive MRI of human gliomas at 3 T using multiecho amine proton chemical exchange saturation transfer spinandgradient echo echoplanar imaging (CESTSAGEEPI). Magn Reson Med. 2018. [Epub ahead of print]

Zhou IY, Wang EF, Cheung JS, Zhang XA, Fulci G, Sun PZ. Quantitative chemical exchange saturation transfer (CEST) MRI of glioma using Image Downsampling Expedited Adaptive Leastsquares (IDEAL) fitting. Sci Rep. 2017;7:84.

Ellingson BM, Bendszus M, Boxerman J, Barboriak D, Erickson BJ, Smits M, Nelson SJ, Gerstner E, Alexander B, Goldmacher G, Wick W, Vogelbaum M, Weller M, Galanis E, KalpathyCramer J, Shankar L, Jacobs P, Pope WB, Yang D, Chung C, Knopp MV, Cha S, van den Bent MJ, Chang S, Yung WK1, Cloughesy TF1, Wen PY1, Gilbert MR1; Jumpstarting Brain Tumor Drug Development Coalition Imaging Standardization Steering Committee. Consensus recommendations for a standardized Brain Tumor Imaging Protocol in clinical trials. Neuro Oncol. 2015;17:1188–1198.

Jenkinson M, Bannister P, Brady M, Smith S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage. 2002;17:825–841.
Research Articles
Download PDF (6.31 MB)
TOMOGRAPHY, September 2018, Volume 4, Issue 3:123137
DOI: 10.18383/j.tom.2018.00017
Improving B0 Correction for pHWeighted Amine Proton Chemical Exchange Saturation Transfer (CEST) Imaging by Use of kMeans Clustering and Lorentzian Estimation
Jingwen Yao^{1}, Dan Ruan^{4}, Catalina Raymond^{1}, Linda M. Liau^{6}, Noriko Salamon^{2}, Whitney B. Pope^{2}, Phioanh L. Nghiemphu^{5}, Albert Lai^{5}, Timothy F. Cloughesy^{5}, Benjamin M. Ellingson^{1}
Abstract
Amine chemical exchange saturation transfer (CEST) echoplanar imaging (EPI) provides unique pH and amino acid MRI contrast, enabling sensitive detection of altered microenvironment properties in various diseases. However, CEST contrast is sensitive to static magnetic field (B_{0}) inhomogeneities. Here we propose 2 new B_{0} correction algorithms for use in correcting pHweighted amine CEST EPI based on kmeans clustering and Lorentzian fitting of CEST data: the iterative downsampling estimation using Lorentzian fitting and the 2stage Lorentzian estimation with 4D polynomial fitting. Higher quality images of asymmetric magnetization transfer ratio (MTR_{asym}) at 3.0 ppm could be obtained with the proposed algorithms than with the existing B_{0} correction methods. In particular, the proposed methods are shown to improve the intertissue consistency, interpatient consistency, and tumor region signaltonoise ratio of MTR_{asym} at 3.0 ppm images, with nonexcessive computation time.
Introduction
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 (MTR_{asym}) = S( − Δω)/S_{0} − S(Δω)/S_{0}, where Δω is the frequency shift between the saturation frequency of interest and the frequency of water proton resonance (or the Zspectrum offset frequency), S(Δω) is the water signal with the application of offresonance saturation pulses, and S_{0} is the reference water signal without the application of saturation pulses. This measure, MTR_{asym}, reduces the direct water saturation and the symmetric macromolecular magnetic transfer effects. However, because calculation of MTR_{asym} requires prior knowledge of water proton resonance frequency, which is directly proportional to the static magnetic field B_{0}, traditional measurements of CEST contrast can be highly sensitive to static magnetic field inhomogeneities (Figure 1).
Figure 1.
Example of asymmetric magnetization transfer ratio (MTR_{asym} at 3.0 ppm) map with and without B_{0} correction. (A) and (C) demonstrate MTR_{asym} (at 3.0 ppm) of the same image slice with and without B_{0} correction, respectively. (B) is the ΔB_{0} map used for correction. Note the higher extent of MTR_{asym} contrast in (A) compared with that after B_{0} correction in (C).
Various approaches have been used to correct B_{0} inhomogeneities. For example, one method involves acquiring additional Zspectral 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 Zspectrum is adjusted accordingly (Figure 2). Although this method is fast and simple, it has several disadvantages including errors introduced because of discretized estimates of ΔB_{0}, restrictions on the detectable range of B_{0} inhomogeneities, and ΔB_{0} estimates are easily affected by signal fluctuation because of noise and other contamination. The frequency resolution of the ΔB_{0} estimates could be improved with fitting and interpolation of Zspectra, including highorder polynomial fitting (6, 8), smoothingsplines 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 Zspectra scan with low B1 amplitude to isolate the effect of direct water saturation.
Figure 2.
Example of Zspectrum B_{0} correction with different methods. Central spectral points correction (CSC) methods (A): The minimum point of the original Zspectrum (blue) is found to be at +0.2 ppm. The whole Zspectrum is shifted so that the minimum point occurs at 0 ppm. The ΔB_{0} map generated by this method is shown on the right. Multiechophase method (B): the ΔB_{0} map is generated from the subtraction of phase images at 2 different echo times, TE1 and TE2. Highorder polynomial (C), smoothingspline (D), singlepool Lorenztian fitting (LF) method (E): the Zspectrum data points are fitted with a ninthorder polynomial function (C)/smoothingspline function (D)/singlepool Lorenztian function (E). The minimum points of the fitting curves are considered to be the B_{0} offsets.
A second method is to use a more traditional approach to B_{0} field mapping through estimating phase differences from 2 gradientecho acquisitions with different echo times, or by using a multiecho acquisition scheme (12). A ΔB_{0} map can be estimated from the difference in the phase measured between 2 different gradientecho times: ΔB_{0} = (ϕ_{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 ΔB_{0} value is known; restrictions on the detectable range of B_{0} 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 B_{0} fluctuations (13).
Figure 3.
Flowchart of the proposed ΔB_{0} estimation methods. The CEST data were preprocessed with motion correction and normalization to reference scans. The normalized CEST data were then used to perform ΔB_{0} estimation with the 2 newly proposed methods. The details of the methods are described in the Methodology.
In this study we introduce and demonstrate the benefits of 2 new methods for B_{0} correction, specifically designed for fast clinical pHweighted amine proton CEST echoplanar imaging (7, 14), by using kmeans clustering combined with Lorentzian fitting (LF) to obtain betterquality images of MTR_{asym} 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 signaltonoise ratio (SNR) of MTR_{asym} at 3.0 ppm between an iterative downsampling estimation (IDE) (15) and a 4D polynomial fitting technique compared with existing approaches.
Methodology
MRI Acquisition
CEST images for all experiments were acquired using a spinandgradient 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 fastexchanging amine protons (7), consisting of three 100millisecond Gaussian pulses with peak amplitude B_{1} = 6 μT were used to sample a total of 29 Zspectral 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 S_{0} 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 wholebrain CESTSAGEEPI and S_{0} images was 7 minutes and 30 seconds. In addition to CESTSAGEEPI acquisition before contrast administration, all patients received anatomic imaging including T2weighted fluidattenuated inversion recovery (FLAIR) images, isotropic 3D T1weighted scans before and after injection of 0.01 mg/kg GdDTPA according to the international standardized brain tumor imaging protocol (16).
Patients
A total of 27 adult primary brain tumor patients (33 scans) with highquality amine CESTSAGEEPI 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 IRBapproved 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.
Table 1.
Patient Characteristics
i] Abbreviation: WHO, World Health Organization.
ΔB_{0} 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 S_{0} reference scan data, resulting in Zspectra data. The Zspectral 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 Zspectral points around the theorized water resonance frequency were compared. The spectral point corresponding to the minimum signal intensity was considered as the B_{0} offset, and the entire Zspectrum was adjusted accordingly (Figure 2A).
ZSpectrum Fitting Correction Method.
Similar to the central spectrum point correction (CSC) method, raw CEST data were first normalized by the S_{0} reference scan, resulting in Zspectra data. Then, for every voxel, the Zspectrum was fitted using a ninthorder polynomial, a smoothing spline, or a singlepool 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 singlepool LF was applied with a constraint of −1 ppm < ΔB_{0} < +1 ppm. After interpolation, the B_{0} offset was estimated as the spectral point corresponding to the minimum signal intensity. The Zspectrum was then adjusted accordingly (Figure 2, C–E).
MultiechoPhase ΔB_{0} Mapping.
B_{0} correction was performed by creating a voxelwise ΔB_{0} map from phase data obtained from the 2 gradient echoes in the S_{0} images, using the following formula: ΔB_{0} = (phase(TE_{2}) − phase(TE_{1})) / γ(TE_{2} − TE_{1}), where TE_{1} and TE_{2} are the echo times of the 2 multiecho images and γ is the gyromagnetic ratio (Figure 2B).
IDE Using LF.
A slicebyslice ΔB_{0} map estimation was achieved iteratively m times through the following steps (Figure 3A).
Step 1.
Raw CEST data at each offset frequency were normalized with respect to the S_{0} reference scan data, resulting in Zspectral measurements.
Step 2.
kmeans clustering (kmeans function, MATLAB release 2017b) was used to cluster the pixels with similar patterns in the Zspectrum. No spatial continuity restriction was applied on the clustering. The number of clusters for each iteration was set to 3^{m}, where m is the iteration number. The original resolution of images can be reached within 7 iterations.
Step 3.
The average Zspectra data were averaged across all pixels within the same cluster.
Step 4.
Mean Zspectra data were fit to the Lorentzian function: c − c / [(x − a)^{2} + b], where a indicates the shift in the Zspectrum center point, b relates to the width of the Lorentzian shape, c and b together determine the yscaling of the function. The fitting started with loosely constrained bounds and narrowed the constraint boundaries after each iteration, m. The constraint bound for ΔB_{0} estimation was set to be ±1/m ppm. If the goodness of fit was low, R^{2} < 0.7, then ΔB_{0} estimation was performed with the conventional CSC method to avoid unrealistic fitting to poor data.
Step 5.
The B_{0} inhomogeneity (ΔB_{0}) estimation (value of a in the Lorentzian fit) was stored for all pixels, then the Zspectra data were updated using the newly estimated ΔB_{0} map.
Step 6.
Determine whether to terminate the iterative process by determining whether the mean square ΔB_{0} estimates for the current iteration was <10% of the sum of the mean square ΔB_{0} estimates over all past iterations, indicating little improvement in estimation by the current iteration. (In practice, the iterative ΔB_{0} estimation typically terminates after 4 or 5 iterations). If this stopping criterion was met, the sum of ΔB_{0} maps from all iterations was determined to be the final ΔB_{0} map. If this criterion was not met, steps 2 through 6 were repeated while the number of kmean clusters were increased and the constraint bounds for ΔB_{0} estimation were tightened.
Step 7.
Lastly, a 3D Gaussian filter (standard deviation = 1) was applied to the final ΔB_{0} map to remove potential erroneous estimations.
TwoStage Lorentzian Estimation With 4D Polynomial Fitting (LFPoly).
A wholebrain 3D ΔB_{0} map estimation was achieved through the following 2stage procedure (Figure 3B).
Stage 1.
Raw CEST data were normalized with respect to the S_{0} reference scan data to calculate the Zspectra. kmeans clustering was used to cluster voxels according to similar patterns in the Zspectrum. The average Zspectrum was calculated for all voxels within the same cluster. Lastly, the mean Zspectrum data were fit with Lorentzian function as described in the previous section.
4D Polynomial Fitting.
ΔB_{0} estimation from stage 1 was fit to a spatial 4D polynomial function (x, y, z spatial coordinates plus the Zspectral dimension). This approach was based on common knowledge that ΔB_{0} is not discontinuous and thus is spatially smooth across the space.
Stage 2.
The Zspectral data were updated using the new estimated ΔB_{0} map following 4D polynomial fitting. Stage 1 was repeated to correct for the residue of the ΔB_{0} map. Lastly, the final ΔB_{0} map estimations was the sum of ΔB_{0} 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 kmeans clusters and the order of 4D polynomial fitting in the LFPoly method, combinations of 8 cluster numbers (10, 20, 50, 100, 200, 300, 400, 500) and 8 polynomial orders (5–12) were used to generate ΔB_{0} estimations and resulting MTR_{asym} 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.
Methods' Comparison
All methods were compared with the following 2 image quality metrics of ΔB_{0} map and resulting MTR_{asym} map (at 3.0 ppm): peak signaltonoise ratio (PSNR) and structural similarity index (SSIM). PSNR was calculated as follows: PSNR = 10log_{10}(peakval^{2}/MSE), where peakval is set to 1 for ΔB_{0} map and 10 for the MTR_{asym} map, and MSE is the mean square error with the reference image. The ΔB_{0} maps and resulting MTR_{asym} images generated by the multiechophase 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; voxelwise Zspectrum fitting using ninthorder polynomial, smoothing spline, and singlepool Lorentzian line fit; IDELF; and 2stage LFPoly.
Medians and standard deviations of the wholebrain MTR_{asym} (at 3.0 ppm) histogram and normalappearing white matter (NAWM) MTR_{asym} (at 3.0 ppm) of all 27 patients (33 scans) were compared among CSC method, IDELF method, multiechophase method, and LFPoly method, by means of Student t test. The MTR_{asym} (at 3.0 ppm) statistics in FLAIR hyperintensity regions were also compared among the aforementioned methods. The SNR of MTR_{asym} in the FLAIR hyperintensity regions was calculated as the ratio between the median MTR_{asym} in FLAIR lesion and the standard deviation of MTR_{asym} 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).
Results
An example of the IDELF ΔB_{0} map estimation is shown in Figure 4. This figure illustrates the ΔB_{0} map estimation from each iteration (Figure 4A) and the sum of ΔB_{0} map estimation after each iteration (Figure 4B). Note that the number of kmeans clusters increased with iteration, and as such, a finer division of brain tissue according to Zspectra data can be observed. In addition, note that the gradient of ΔB_{0} also becomes smoother in later iterations. An example of mean square (MS) ΔB_{0} map – iteration number relationship from one slice is plotted in Figure 4. The MS value of ΔB_{0} map from each iteration decreases with increase in the iteration number, indicating less extra information about ΔB_{0} map from later iterations. The stopping criterion of MS(ΔB_{0}) for the current iteration being smaller than 10% of the sum of MS(ΔB_{0}) 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).
Figure 4.
Example of Iterative downsampling estimation (IDELF) of ΔB_{0} map. (A) shows the ΔB_{0} map estimation from iterations 1–7 with kmeans clustering of increasing cluster numbers. The mean square (MS) of the ΔB_{0} map estimation decreases with the increase in iteration number, indicating less extra information of the ΔB_{0} map obtained from each iteration. An example of the relationship between MS values and iteration numbers for one slice of the brain is plotted. (B) demonstrates the sum of ΔB_{0} estimations after iterations 1–7. The ΔB_{0} estimations become smoother and contain more details as the iteration number increases. The last image shows the final ΔB_{0} map after Gaussian filtering.
The number of kmeans 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 ΔB_{0} estimation. The PSNR and SSIM metrics of ΔB_{0} maps and MTR_{asym} 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 kmeans clusters and the increasing polynomial fitting order. However, the effects of kmeans cluster number appeared more influential than those of the polynomial fitting order, and in most circumstances, this increase plateaued when the number of kmeans 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 kmeans clusters and eighthorder 4D polynomial fitting and 500 kmeans clusters and 12thorder 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).
Figure 5.
Image quality metrics and computation time for the 2stage Lorentzian estimation with 4D polynomial fitting (LFPoly) method with differing numbers of kmeans clusters (range from 10 to 500) and orders of 4D polynomial fitting (range from 5 to 12). The figure shows the PSNR (A) and SSIM (B) metrics of ΔB_{0} maps and MTR_{asym} images (C, D) generated with different combinations of parameters. All image quality metrics are normalized to the highest values of all combinations of parameters. The average of the 4 metrics is plotted in (E). The computation time required for each combination of parameters is plotted in (F). The blue and red circles correspond to the cases with 100 kmeans clusters + eighthorder 4D polynomial fitting and 500 kmeans clusters + 12thorder 4D polynomial fitting, respectively. The 2 parameter combinations were used for further comparison with other methods.
An example slice of estimated ΔB_{0} map and resulting MTR_{asym} image generated by different correction methods is illustrated in Figure 6. The quantitative assessments of ΔB_{0} estimation methods for 7 patients are plotted in Figure 7, including image quality metrics PSNR and SSIM of ΔB_{0} maps and MTR_{asym} images, computation time, and standard deviation of NAWM MTR_{asym}. The image quality metrics and the standard deviation of NAWM MTR_{asym} show consistent results. Among all methods except for multiechophase method, the proposed IDELF method and the LFPoly method perform the best in terms of having the highest PSNR and SSIM values and lowest NAWM MTR_{asym} standard deviations, while the ninthorder polynomial and the smoothingspline fittings of Zspectra have the worst performance. Consistent with the quantitative metrics, the proposed IDELF method and the LFPoly method generated MTR_{asym} images that are visually less noisy in Figure 6, while the ninthorder polynomial and the smoothingspline fitting methods generated highly noisy MTR_{asym} images, particularly around the nasal sinuses, indicating very inaccurate estimation of ΔB_{0}.
Figure 6.
ΔB_{0} maps (A) and resulting MTR_{asym} images (B) generated by different methods. IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting. Black arrow: MTR_{asym} (at 3.0 ppm) contrast in tumor area; white arrow: region with severe B_{0} inhomogeneity.
Figure 7.
Image quality metrics for ΔB_{0} maps and resulting MTR_{asym} (at 3.0 ppm) images generated by different methods (data from 7 patients). The PSNR and SSIM metrics of ΔB_{0} maps and MTR_{asym} (at 3.0 ppm) images are shown in (A–D). The computation times of different methods are plotted in (E). (F) shows the standard deviation of normal appearing white matter (NAWM) MTR_{asym} (at 3.0 ppm). Subplots were included in (E) and (F) to zoom in the metric scale for selected methods. Triangle markers of different colors represent data points from different patients. CSC: central spectrum points correction method; IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 kmeans clusters and eighthorder polynomial and (b) parameters of 500 kmeans clusters and 12thorder polynomial. Box plot: the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles.
The computation time for wholebrain ΔB_{0} estimation using different methods varies with orders of magnitude. The CSC method and the multiechophase method have near realtime execution time. The LFPoly method has the computation time of ∼≤1 minute. The IDELF method, as well as the ninthorder polynomial and the smoothingspline fitting methods, requires a computation time of 10–15 minutes. Voxelbyvoxel singlepool Lorentzian line fitting requires the longest computation time of ∼40 minutes to 1 hour per patient.
The CSC method, the multiechophase method, and the proposed IDELF and LFPoly 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 IDELF method and the LFPoly method show significantly higher intertissue consistency (lower standard deviation of NAWM MTR_{asym} at 3.0 ppm) than the conventional CSC method and the multiechophase method (IDELF: P = 2.03 × 10^{−18} compared with the CSC method and P = 1.10 × 10^{−4} compared with the the multiechophase method; LFPoly(a): P = 9.02 × 10^{−19} compared with the CSC method and P = 3.03 × 10^{−4} compared with the multiechophase method). There is no significant difference in the standard deviation of NAWM MTR_{asym} at 3.0 ppm between IDELF method and the LFPoly methods (P = .402/.939). The proposed methods also show higher interpatient consistency (lower standard deviation of median NAWM MTR_{asym} at 3.0 ppm) than the CSC and the multiechophase methods. The standard deviation of median NAWM MTR_{asym} is 0.2953% for the IDELF method and 0.2812%/0.2666% for the LFPoly methods, compared with 0.3025% for the CSC method, and 2.0011% for the multiechophase method. Between the LFPoly methods with different parameter choices, a slightly lower standard deviation of NAWM MTR_{asym} and a slightly lower standard deviation of median NAWM MTR_{asym} can be observed in the case of 500 kmeans clusters and 12thorder polynomial fitting compared with 100 kmeans clusters and eighthorder polynomial fitting, but the difference is not significant (P = .2390). A more detailed comparison of NAWM MTR_{asym} statistics can be found in Table 2.
Figure 8.
MTR_{asym} statistics with selected B_{0} correction methods (data from 27 patients). Standard deviations (A) and medians of NAWM MTR_{asym} (at 3.0 ppm) (B) are plotted. The IDELF method shows the lowest standard deviation of NAWM MTR_{asym} (at 3.0 ppm) and the lowest intersubject deviation of median NAWM MTR_{asym} (at 3.0 ppm) (***P < .001). The computation time of wholebrain ΔB_{0} estimation with different methods is shown in (C). The CSC method and the multiechophase method require the least computation time (median computation time, 1.0 and 0.7 seconds, respectively), while the IDELF method requires the most (median computation time, 697.2 seconds). Examples of wholebrain MTR_{asym} (at 3.0 ppm) histograms of 4 patients' data is demonstrated in (D). Generally speaking, the IDELF method and the LFPoly method lead to narrower distribution, while multiechophase method shifts the distribution in some of the cases. CSC: central spectrum points correction method; IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 kmeans clusters and 8th order polynomial, (b) with parameters of 500 kmeans clusters and 12th order polynomial. Box plot: the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles.
Figure 9.
Comparison of ROI dependence on MTRasym statistics. FLAIR hyperintensity region of interest (ROI) MTR_{asym} (at 3.0 ppm) SNR (calculated as median of MTR_{asym} / standard deviation of NAWM MTR_{asym}) in presurgery scans of patients with glioma is shown in (A). The IDELF method and the LFPoly method shows significantly higher SNR than the conventional CSC method (***P < .001). The phase method also shows higher SNR. However, the intersubject deviation is large, similar to the case in NAWM MTR_{asym}. Comparison of MTR_{asym} (at 3.0 ppm) median values in NAWM, FLAIR lesion in presurgery scans of patients with glioma, FLAIR lesion in postsurgery scans (postOP) of patients with glioma, with IDELF B_{0} correction method is shown in (B). NAWM regions have significantly lower MTR_{asym} median compared with presurgery glioma FLAIR lesions (***P < .001) and with postsurgery FLAIR lesions (***P < .001). Postsurgery FLAIR lesions have significantly lower MTR_{asym} median than presurgery glioma FLAIR lesions (**P < .01). CSC: central spectrum points correction method; IDELF: iterative downsampling (kmeans clustering) estimation with Lorentzian fitting; LFPoly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 kmeans clusters and 8^{th} order polynomial, (b) with parameters of 500 kmeans clusters and 12^{th} order polynomial. Box plot: the central mark indicates the median, and the bottom and top edges of the box indicate the 25^{th} and 75^{th} percentiles.
Table 2.
NAWM MTR_{asym} (at 3.0 ppm) Statistics
i] ***P < .001.
In Figure 9, A and B, the FLAIR hyperintensity lesion MTR_{asym} (at 3.0 ppm) with different B_{0} correction methods are shown. The IDELF method and the LFPoly methods show significantly higher FLAIR lesion MTR_{asym} SNR than the conventional CSC method (IDELF: P = 2.74 × 10^{−4} compared with the CSC method; LFkPoly(a): P = 2.08 × 10^{−5} compared with the CSC method; LFkPoly(b): P = 3.04 × 10^{−5} compared with the CSC method). Although the multiechophase 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 MTR_{asym} statistics can be found in Table 3. With the proposed IDELF method, a significant difference in median MTR_{asym} values can be observed in different ROIs and different groups of patients. Presurgery glioma FLAIR lesions (12 scans) have significantly higher MTR_{asym} 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 MTR_{asym} median than presurgery glioma FLAIR lesions (P = .0022).
Table 3.
Presurgery FLAIR Lesion MTR_{asym} (at 3.0 ppm) Statistics in Patients With Glioma
i] ***P < .001.
The computation time of the selected methods (Figure 8C) was consistent with the preliminary results from 7 patient data. The IDELF method requires the longest computation time (680.4 seconds ± 342.2 seconds per patient), while the LFPoly method can provide wholebrain ΔB_{0} estimation in ∼1 min, 28.6 seconds ± 3.3 seconds for LFPoly(a) and 65.9 seconds ± 5.3 seconds for LFPoly(b). The CSC method and the multiechophase method have near realtime computation time.
In Figure 8D, the histograms of wholebrain MTR_{asym} using different B_{0} correction methods were plotted for 4 patients as examples. Compared with the CSC method, the multiechophase method, the IDELF method, and the LFPoly method all generated narrower distribution of MTR_{asym} (at 3.0 ppm), possibly indicating improved intertissue consistency. However, the multiechophase method shifted the distribution noticeably in some of the patients. The proposed IDELF and LFPoly methods shifted the distribution slightly toward the negative direction, which could be explained by the better ΔB_{0} estimation close to the sinus regions, alleviating the overestimation of MTR_{asym} (at 3.0 ppm) in areas with very high ΔB_{0} (Figure 6). No difference of wholebrain MTR_{asym} distribution has been observed between the IDELF and LFPoly methods.
An example of wholebrain ΔB_{0} correction with multiechophase method and IDELF method can be seen in Figure 10 to better understand the deviation of ΔB_{0} estimation between the 2 methods. Several differences have been observed:
The proposed IDELF method avoided overestimation of ΔB_{0} in slice 16 and slice 22, indicated by the more consistent white matter MTR_{asym} across the slices.
The IDELF method does not suffer from phaseunwrapping errors, which can be observed in slice 4 for the multiechophase method, indicated by the abnormal low MTR_{asym} signal in the right temporal lobe.
Figure 10.
Wholebrain ΔB_{0} maps and resulting MTR_{asym} (at 3.0 ppm) maps generated by the multiechophase method (A, B) and the IDELF method (C, D). The MTR_{asym} (at 3.0 ppm) maps generated by the IDELF method show more consistency across the slices compared with the multiechophase method.
Discussion
This study explored the application of 2 newly proposed methods using kmeans clustering combined with LF in B_{0} 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 ΔB_{0}. The results clearly show that the proposed methods outperform the existing methods (CSC method, polynomial/smoothingspline/singlepool LF, multiechophase method) within the context of clinical CEST data that are often undersampled along the Zspectrum. The proposed methods can serve as promising new approaches for improving CEST image B_{0} 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 Zspectra spectral points compared with research scans, where more number of averages and spectral points can be acquired. Even with accelerated CEST sequences, such as CESTEPISAGE (14), only one measurement of 29 Zspectral points was possible within clinically allowed time (<10 minute) for a wholebraincoverage CEST imaging. This limitation restricted the use of many B_{0} 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 highorder polynomial fitting and smoothingspline fitting of Zspectra. The poor performance shown in Figure 7 was likely because of the error induced in fitting noisy data and unevenly sampled Zspectral points (Figure 2). The uneven sampling of the Zspectrum particularly impairs the fitting accuracy of the polynomial method when B_{0} inhomogeneity is large (outside the range of central spectral points). This can be observed in Figure 6. Although voxelbyvoxel singlepool LF performs better, as it uses the global information of Zspectra; this method requires unrealistic computation time for wholebrain analysis (40 minutes to 1 hour per patient) and may still suffer from corrupt Zspectral data.
Compared with the conventional CSC method and the multiechophase method, the 2 newly proposed approaches have the advantages of both methods, while avoiding many of their drawbacks. For example, the proposed IDELF and LFPoly methods use the acquired CEST spectra data to perform inline correction, avoiding the potential issue with ΔB_{0} map fluctuations and added acquisition time. These methods also avoid potential phaseunwrapping 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 B_{0} 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 ΔB_{0} map through an iterative downsampling method in IDELF or by 4D spatial polynomial fitting of ΔB_{0} map in LFPoly, avoiding potential erroneous values that may arise from other methods.
The improved performance of ΔB_{0} estimation using kmeans clustering combined with LF could be explained by the “denoising” effect on the resulting Zspectrum data. kmeans Clustering of 4D data cluster voxels with similar patterns and features in the Zspectral data and thus often group similar tissue types together. The averaging of Zspectral data within the same cluster greatly improves the SNR of the resulting Zspectral data, thus allowing better estimation of ΔB_{0} and faster convergence of LF.
In the proposed methods, a single Lorentzian line shape was used to fit the underlying Zspectral data. This may have caused potential issues with overestimating ΔB_{0} in voxels with lower values of positive Zspectrum points (ie, higher MTR_{asym}). The weighting of central spectrum points (from −0.3 ppm to +0.3 ppm) was increased to mitigate this issue, but no overestimation of ΔB_{0} 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 Zspectra, particularly when the asymmetry approaches free water proton resonance, for example, when investigating the hydroxyl (OH) proton pool.
The two proposed methods for ΔB_{0} estimation, IDELF and LFPoly, show similar performance resulting in improved image quality in MTR_{asym}, reduced intrasubject and intersubject NAWM MTR_{asym} deviations, and improved SNR in FLAIR lesion MTR_{asym}. The LFPoly accounts for the spatial smoothness of ΔB_{0} inhomogeneities with the 4D polynomial fitting, which at the same time could cause less accurate estimations of local ΔB_{0} details. However, the much shorter computation time of LFPoly (∼30 seconds) compared with IDELF (∼10 minutes) makes it a better choice for fast generation of MTR_{asym} images in clinical use, while IDELF may be considered for a more indepth quantitative analysis of CEST data when long reconstruction and postprocessing times are not an issue.
No significant difference was observed between the performance of LFPoly methods with the 2 different parameter settings (100 kmeans clusters + eighthorder polynomial fitting and 500 kmeans clusters + 12^{th}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 Zspectral 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 kmeans clusters and eighthorder 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 Zspectral 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 highquality ΔB_{0} 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 ΔB_{0} map.
In conclusion, we proposed 2 new methods for CEST B_{0} correction using kmeans clustering combined with LF for improving B_{0} correction for fast clinical pHweighted CEST imaging: the iterative downsampling estimation using Lorentzian fitting (IDELF) and the 2stage Lorentzian estimation with 4D polynomial fitting (LFPoly). The proposed methods were shown to improve the image quality, intertissue consistency, interpatient consistency, and lesion SNR of MTR_{asym} (at 3.0 ppm) images.
Notes
[4] Abbreviations:
CEST
Chemical exchange saturation transfer
EPI
echoplanar imaging
MTR_{asym}
asymmetric magnetization transfer ratio
MRI
magnetic resonance imaging
LF
Lorentzian fitting
MR
magnetic resonance
IDE
iterative downsampling estimation
SNR
signaltonoise ratio
SAGE
spinandgradient echo
FLAIR
fluidattenuated inversion recovery
PSNR
peak signaltonoise ratio
SSIM
structural similarity index
ROIs
regions of interest
MS
mean square
WASSR
water saturation shift referencing
Acknowledgements
We gratefully acknowledge the support of the American Cancer Society (ACS), Art of the Brain, University of California Research Coordinating Committee, UCLA Jonsson Comprehensive Cancer Center Seed Grant, UCLA SPORE, and NIH/NCI Funding. The study was funded by the American Cancer Society (ACS) Research Scholar Grant (RSG1500301CCE) (Ellingson); Art of the Brain (Cloughesy); University of California Research Coordinating Committee (Ellingson); UCLA Jonsson Comprehensive Cancer Center Seed Grant (Ellingson); UCLA SPORE in Brain Cancer (NIH/NCI 1P50CA21101501A1) (Ellingson, Liau, Nghiemphu, Lai, Pope, Cloughesy); and NIH/NCI 1R21CA22375701 (Ellingson).
Disclosures: No disclosures to report.
Conflict of Interest: The authors have no conflict of interest to declare.
References
Journal Information
Journal ID (nlmta): tom
Journal ID (publisherid): TOMOG
Title: Tomography
Subtitle: A Journal for Imaging Research
Abbreviated Title: Tomog.
ISSN (print): 23791381
ISSN (electronic): 2379139X
Publisher: Grapho Publications, LLC (Ann Abor, Michigan)
Article Information
Self URI: tom00318000123.pdf
Copyright statement: © 2018 The Authors. Published by Grapho Publications, LLC
Copyright: 2018, Grapho Publications, LLC
License (openaccess, http://creativecommons.org/licenses/byncnd/4.0/):
This is an open access article under the CC BYNCND license (http://creativecommons.org/licenses/byncnd/4.0/).
Publication date (print): September 2018
Volume: 4
Issue: 3
Pages: 123137
Publisher ID: TOMO201800017
DOI: 10.18383/j.tom.2018.00017
PDF
Download the article PDF (6.31 MB)
Download the full issue PDF (9.12 MB)
Mobileready Flipbook
View the full issue as a flipbook (Desktop and Mobileready)