Research Articles

Download PDF (6.31 MB)

TOMOGRAPHY, September 2018, Volume 4, Issue 3:123-137
DOI: 10.18383/j.tom.2018.00017

Improving B0 Correction for pH-Weighted Amine Proton Chemical Exchange Saturation Transfer (CEST) Imaging by Use of k-Means Clustering and Lorentzian Estimation

Jingwen Yao1, Dan Ruan4, Catalina Raymond1, Linda M. Liau6, Noriko Salamon2, Whitney B. Pope2, Phioanh L. Nghiemphu5, Albert Lai5, Timothy F. Cloughesy5, Benjamin M. Ellingson1

1UCLA Brain Tumor Imaging Laboratory (BTIL), Center for Computer Vision and Imaging Biomarkers, David Geffen School of Medicine, University of California Los Angeles, Los Angeles, CA;2Department of Radiological Sciences, David Geffen School of Medicine, University of California Los Angeles, Los Angeles, CA;3Department of Bioengineering, Henry Samueli School of Engineering and Applied Science, University of California Los Angeles, Los Angeles, CA;Departments of 4Radiation Oncology,5Neurology, and6Neurosurgery, David Geffen School of Medicine, University of California Los Angeles, Los Angeles, CA

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 (B0) inhomogeneities. Here we propose 2 new B0 correction algorithms for use in correcting pH-weighted amine CEST EPI based on k-means clustering and Lorentzian fitting of CEST data: the iterative downsampling estimation using Lorentzian fitting and the 2-stage Lorentzian estimation with 4D polynomial fitting. Higher quality images of asymmetric magnetization transfer ratio (MTRasym) at 3.0 ppm could be obtained with the proposed algorithms than with the existing B0 correction methods. In particular, the proposed methods are shown to improve the intertissue consistency, interpatient consistency, and tumor region signal-to-noise ratio of MTRasym 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 (MTRasym) = S( − Δω)/S0S(Δω)/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).

Figure 1.

Example of asymmetric magnetization transfer ratio (MTRasym at 3.0 ppm) map with and without B0 correction. (A) and (C) demonstrate MTRasym (at 3.0 ppm) of the same image slice with and without B0 correction, respectively. (B) is the ΔB0 map used for correction. Note the higher extent of MTRasym contrast in (A) compared with that after B0 correction in (C).

media/vol4/issue3/images/tom0031801200001.jpg

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.

Figure 2.

Example of Z-spectrum B0 correction with different methods. Central spectral points correction (CSC) methods (A): The minimum point of the original Z-spectrum (blue) is found to be at +0.2 ppm. The whole Z-spectrum is shifted so that the minimum point occurs at 0 ppm. The ΔB0 map generated by this method is shown on the right. Multiecho-phase method (B): the ΔB0 map is generated from the subtraction of phase images at 2 different echo times, TE1 and TE2. High-order polynomial (C), smoothing-spline (D), single-pool Lorenztian fitting (LF) method (E): the Z-spectrum data points are fitted with a ninth-order polynomial function (C)/smoothing-spline function (D)/single-pool Lorenztian function (E). The minimum points of the fitting curves are considered to be the B0 offsets.

media/vol4/issue3/images/tom0031801200002.jpg

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).

Figure 3.

Flowchart of the proposed ΔB0 estimation methods. The CEST data were preprocessed with motion correction and normalization to reference scans. The normalized CEST data were then used to perform ΔB0 estimation with the 2 newly proposed methods. The details of the methods are described in the Methodology.

media/vol4/issue3/images/tom0031801200003.jpg

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.

Methodology

MRI Acquisition

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).

Patients

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.

Table 1.

Patient Characteristics

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%)

i] Abbreviation: WHO, World Health Organization.

Δ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).

Step 1.

Raw CEST data at each offset frequency were normalized with respect to the S0 reference scan data, resulting in Z-spectral measurements.

Step 2.

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.

Step 3.

The average Z-spectra data were averaged across all pixels within the same cluster.

Step 4.

Mean Z-spectra data were fit to the Lorentzian function: cc / [(xa)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.

Step 5.

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.

Step 6.

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.

Step 7.

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).

Stage 1.

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.

Stage 2.

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.

Methods' Comparison

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).

Results

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).

Figure 4.

Example of Iterative down-sampling estimation (IDE-LF) of ΔB0 map. (A) shows the ΔB0 map estimation from iterations 1–7 with k-means clustering of increasing cluster numbers. The mean square (MS) of the ΔB0 map estimation decreases with the increase in iteration number, indicating less extra information of the ΔB0 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 ΔB0 estimations after iterations 1–7. The ΔB0 estimations become smoother and contain more details as the iteration number increases. The last image shows the final ΔB0 map after Gaussian filtering.

media/vol4/issue3/images/tom0031801200004.jpg

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).

Figure 5.

Image quality metrics and computation time for the 2-stage Lorentzian estimation with 4D polynomial fitting (LF-Poly) method with differing numbers of k-means 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 ΔB0 maps and MTRasym 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 k-means clusters + eighth-order 4D polynomial fitting and 500 k-means clusters + 12th-order 4D polynomial fitting, respectively. The 2 parameter combinations were used for further comparison with other methods.

media/vol4/issue3/images/tom0031801200005.jpg

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.

Figure 6.

ΔB0 maps (A) and resulting MTRasym images (B) generated by different methods. IDE-LF: iterative down-sampling (k-means clustering) estimation with Lorentzian fitting; LF-Poly: two stage Lorentzian estimation with 4D polynomial fitting. Black arrow: MTRasym (at 3.0 ppm) contrast in tumor area; white arrow: region with severe B0 inhomogeneity.

media/vol4/issue3/images/tom0031801200006.jpg
Figure 7.

Image quality metrics for ΔB0 maps and resulting MTRasym (at 3.0 ppm) images generated by different methods (data from 7 patients). The PSNR and SSIM metrics of ΔB0 maps and MTRasym (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) MTRasym (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; IDE-LF: iterative down-sampling (k-means clustering) estimation with Lorentzian fitting; LF-Poly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 k-means clusters and eighth-order polynomial and (b) parameters of 500 k-means 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.

media/vol4/issue3/images/tom0031801200007.jpg

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.

Figure 8.

MTRasym statistics with selected B0 correction methods (data from 27 patients). Standard deviations (A) and medians of NAWM MTRasym (at 3.0 ppm) (B) are plotted. The IDE-LF method shows the lowest standard deviation of NAWM MTRasym (at 3.0 ppm) and the lowest intersubject deviation of median NAWM MTRasym (at 3.0 ppm) (***P < .001). The computation time of whole-brain ΔB0 estimation with different methods is shown in (C). The CSC method and the multiecho-phase method require the least computation time (median computation time, 1.0 and 0.7 seconds, respectively), while the IDE-LF method requires the most (median computation time, 697.2 seconds). Examples of whole-brain MTRasym (at 3.0 ppm) histograms of 4 patients' data is demonstrated in (D). Generally speaking, the IDE-LF method and the LF-Poly method lead to narrower distribution, while multiecho-phase method shifts the distribution in some of the cases. CSC: central spectrum points correction method; IDE-LF: iterative down-sampling (k-means clustering) estimation with Lorentzian fitting; LF-Poly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 k-means clusters and 8th order polynomial, (b) with parameters of 500 k-means 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.

media/vol4/issue3/images/tom0031801200008.jpg
Figure 9.

Comparison of ROI dependence on MTRasym statistics. FLAIR hyperintensity region of interest (ROI) MTRasym (at 3.0 ppm) SNR (calculated as median of MTRasym / standard deviation of NAWM MTRasym) in presurgery scans of patients with glioma is shown in (A). The IDE-LF method and the LF-Poly 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 MTRasym. Comparison of MTRasym (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 IDE-LF B0 correction method is shown in (B). NAWM regions have significantly lower MTRasym median compared with presurgery glioma FLAIR lesions (***P < .001) and with postsurgery FLAIR lesions (***P < .001). Postsurgery FLAIR lesions have significantly lower MTRasym median than presurgery glioma FLAIR lesions (**P < .01). CSC: central spectrum points correction method; IDE-LF: iterative down-sampling (k-means clustering) estimation with Lorentzian fitting; LF-Poly: two stage Lorentzian estimation with 4D polynomial fitting, (a) with parameters of 100 k-means clusters and 8th order polynomial, (b) with parameters of 500 k-means 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.

Table 2.

NAWM MTRasym (at 3.0 ppm) Statistics

Method Comparison(NAWM MTRasym SD) CSCP-value Multi-echo PhaseP-value IDE-LFP-value LF-Poly(a)P-value LF-Poly(b)P-value
CSC N/A
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%)

i] ***P < .001.

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).

Table 3.

Presurgery FLAIR Lesion MTRasym (at 3.0 ppm) Statistics in Patients With Glioma

Method Comparison(FLAIR Lesion SNR) CSCP-value Multiecho PhaseP-value IDE-LFP-value LF-Poly(a)P-value LF-Poly(b)P-value
CSC N/A
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)

i] ***P < .001.

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:

  1. 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.

  2. 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.

Figure 10.

Whole-brain ΔB0 maps and resulting MTRasym (at 3.0 ppm) maps generated by the multiecho-phase method (A, B) and the IDE-LF method (C, D). The MTRasym (at 3.0 ppm) maps generated by the IDE-LF method show more consistency across the slices compared with the multiecho-phase method.

media/vol4/issue3/images/tom0031801200009.jpg

Discussion

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.

Notes

[4] Abbreviations:

CEST

Chemical exchange saturation transfer

EPI

echoplanar imaging

MTRasym

asymmetric magnetization transfer ratio

MRI

magnetic resonance imaging

LF

Lorentzian fitting

MR

magnetic resonance

IDE

iterative downsampling estimation

SNR

signal-to-noise ratio

SAGE

spin-and-gradient echo

FLAIR

fluid-attenuated inversion recovery

PSNR

peak signal-to-noise 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 (RSG-15-003-01-CCE) (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 1P50CA211015-01A1) (Ellingson, Liau, Nghiemphu, Lai, Pope, Cloughesy); and NIH/NCI 1R21CA223757-01 (Ellingson).

Disclosures: No disclosures to report.

Conflict of Interest: The authors have no conflict of interest to declare.

References

  1.  
    Zaiss M, Bachert P. Exchange-dependent relaxation in the rotating frame for slow and intermediate exchange-modeling off-resonant spin-lock and chemical exchange saturation transfer. NMR Biomed. 2013;26:507–518.
  2.  
    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.
  3.  
    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.
  4.  
    Harris RJ, Cloughesy TF, Liau LM, Prins RM, Antonios JP, Li D, Yong WH, Pope WB, Lai A, Nghiemphu PL, Ellingson BM. pH-Weighted molecular imaging of gliomas using amine chemical exchange saturation transfer MRI. Neuro Oncol. 2015;17:1514–1524.
  5.  
    Langereis S, Keupp J, van Velthoven JL, de Roos IH, Burdinski D, Pikkemaat JA, Grüll H. A temperature-sensitive liposomal 1H CEST and 19F contrast agent for MR image-guided drug delivery. J Am Chem Soc. 2009;131:1380–1381.
  6.  
    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.
  7.  
    Harris RJ, Cloughesy TF, Liau LM, Nghiemphu PL, Lai A, Pope WB, Ellingson BM. Simulation, phantom validation, and clinical evaluation of fast pH-weighted molecular imaging using amine chemical exchange saturation transfer echo planar imaging (CEST-EPI) in glioma at 3 T. NMR Biomed. 2016;29:1563–1576.
  8.  
    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.
  9.  
    Stancanello J, Terreno E, Castelli DD, Cabella C, Uggeri F, Aime S. Development and validation of a smoothing-splines-based correction method for improving the analysis of CEST-MR images. Contrast Media Mol Imaging. 2008;3:136–149.
  10.  
    Zaiss M, Schmitt B, Bachert P. Quantitative separation of CEST effect from magnetization transfer and spillover effects by Lorentzian-line-fit analysis of z-spectra. J Magn Reson. 2011;211:149–155.
  11.  
    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.
  12.  
    Sun PZ, Farrar CT, Sorensen AG. Correction for artifacts induced by B-0 and B-1 field inhomogeneities in pH-sensitive chemical exchange saturation transfer (CEST) Imaging. Magn Reson Med. 2007;58:1207–1215.
  13.  
    Durand E, van de Moortele PF, Pachot-Clouard M, Le Bihan D. Artifact due to B-o fluctuations in fMRI: correction using the k-space central line. Magn Reson Med. 2001;46:198–201.
  14.  
    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 pH-sensitive and oxygen-sensitive MRI of human gliomas at 3 T using multi-echo amine proton chemical exchange saturation transfer spin-and-gradient echo echo-planar imaging (CEST-SAGE-EPI). Magn Reson Med. 2018. [Epub ahead of print]
  15.  
    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 Least-squares (IDEAL) fitting. Sci Rep. 2017;7:84.
  16.  
    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, Kalpathy-Cramer 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.
  17.  
    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.

PDF

Download the article PDF (6.31 MB)

Download the full issue PDF (9.12 MB)

Mobile-ready Flipbook

View the full issue as a flipbook (Desktop and Mobile-ready)