Whole brain snapshot CEST at 3T using 3D-EPI: Aiming for speed, volume, and homogeneity

Purpose: CEST MRI enables imaging of distributions of low-concentrated metabolites as well as proteins and peptides and their alterations in diseases. CEST examinations often suffer from low spatial resolution, long acquisition times, and concomitant motion artifacts. This work aims to maximize both resolution and volume coverage with a 3D-EPI snapshot


| INTRODUCTION
CEST MRI allows indirect detection of low concentrated solutes. CEST uses two properties of these low concentrated solutes: their chemical exchange with the bulk water proton pool, and their chemical shift relative to water by a certain frequency offset δω. In contrast to purely anatomical MRI, a frequency selective presaturation module with RF irradiation at δω precedes each image readout. Given presaturation at different frequencies, the resulting bulk water signals may be combined to a Z-spectrum, 1 which contains information on all interacting metabolites resonating in the covered frequency range. To extract all information, the Z-spectrum should be sampled densely, and the presaturation and readout block are thus repeated 10 to 100 times. The major problem of the combination of CEST presaturation and classic line-byline MRI readout is that the repeated excitations will successively destroy the prepared magnetization state. The number of excitations can be reduced by k-space segmentation or by a small matrix size (small FOV or low resolution). Recently, a snapshot CEST approach with centric reordered 3D-gradient echo (GRE) readout scheme was proposed that was optimized for maximized number of k-space lines being acquired following a single presaturation. 2 Still, this GRE-based method was limited to 16 slices and a resolution of 1.7 × 1.7 × 5 mm 3 . When extending the volume coverage to whole brain, even at 3T CEST effects become sensitive to the inhomogeneity of B 1 RF irradiation. To correct for spatial inhomogeneity in B 1 , repeated measurements at different nominal B 1 values are needed. 3 In this work, encouraged by recent findings at 7T, 4 we extended a 3D-EPI readout to allow semielliptical scanning along with CAIPIRINHA 5 undersampling. It was then optimized with centric reordering and nonselective excitation for whole-brain snapshot CEST MRI at 3T. This increased temporal SNR (tSNR) by 20% compared to the 7T protocol applied at 3T. In addition, a novel saturation scheme at low power B 1 was applied to minimize side bands in the Z-spectrum. Together with the proposed postprocessing, this enables multipool CEST MRI with whole brain coverage and isotropic resolution of 1.8 mm in less than 4:30 minutes for a fully sampled Z-spectrum at a clinical MR system.

| MR examinations
All examinations were performed at a clinical whole-body MR system (3T Magnetom Prisma fit , Siemens Healthineers, Germany). The vendor-provided 1Tx/64Rx-channel head/ neck coil was used for MR image acquisition on six healthy subjects with written informed consent and approval by the local ethics committee.
The CEST presaturation module was optimized for relayed nuclear Overhauser enhancement (rNOE) and amide proton transfer (APT) contrast by modifying a scheme suggested by Deshmane et al. 6 To reduce sidebands and increase selectivity, rather long Gaussian pulses of 16 × 100 ms were applied in the present study at a duty cycle of 50% and B 1 = 0.65 µT. The offset distribution can be found in the Supporting Information. To correct for B 1 inhomogeneity, data were acquired at [0.75, 1.00, 1.25] times the nominal B 1 by adjusting the reference voltage. This covered 95% of all measured relative B 1 values in vivo using the body coil as transmit coil. For imaging, both a 3D-GRE 2 and the 3D-EPI 7 sequence were investigated with nonselective excitation applied in both cases.

| MR readout parameters
The 3D-EPI readout was accelerated using parallel acquisition with CAIPIRINHA acceleration 1 × 6 (shift = 2, 36 × 36 reference lines) along both phase-encoding directions (PE, 3D) and 6/8 partial Fourier along the first PE direction. Semielliptical k-space sampling, 7 along with centric reordering of the 3D phase-encoding steps, were used. A non-blipped CAIPIRINHA sampling strategy for high spatiotemporal resolution 8 was chosen for which each CAIPIRINHA offset was realized with a separate excitation instead of 3D gradient blips. This resulted in a maximum EPI-factor of 32 and 45 excitations for the readout of 1201 k-space lines per 3D volume at a nominal matrix size of 144 × 126 × 88 (RO × PE × 3D: head-foot × anteroposterior × left-right). A readout bandwidth of 1930 Hz/pixel with TE/echo train length = 11.0/24.3 ms lead to a total readout duration of 1.2 s for an isotropic nominal resolution of 1.8 mm at a FOV of 256 × 224 × 156 mm 3 . The nominal excitation flip angle (FA) was set to 15 degrees. Nonselective binomial-11 water excitation was applied instead of fat saturation to minimize chemical shift artifacts and maximize readout efficiency. 9 The acquisition duration was 4.3 s per presaturation offset, plus an additional 12 s recovery time for the unsaturated M 0 image. Fifty-seven offsets and one M 0 image could be measured in 4:22 min, including the CAIPIRINHA reference lines.
The 3D-GRE readout was accelerated using GRAPPA 10 3 × 2 in both phase-encoding directions. Further undersampling was achieved by applying partial Fourier (6/8 along both phase-encoding directions) and omitting the corners of k-space (elliptical sampling). The number of acquired k-space lines for a nominal matrix size of 96 × 78 × 72 (RO × PE × 3D: head-foot × anteroposterior × left-right) was therefore 487, and thus only 15% above the maximum of 424 lines derived according to the work of Zaiss et al 2 (T 1,WM = 950 ms 11 ; FA = 5°; TR = 3.1 ms). The whole readout duration was 1.7 s. These optimized settings of the 3D-GRE (TE = 1.3 ms; TR = 3.1 ms; bandwidth = 660 Hz/pixel; FA = 5°) 12 allowed to acquire one presaturation offset in 5.1 s (additional relaxation in case of M 0 image) with a nominal spatial resolution of 2.34 mm isotropic.
Field mapping was performed applying the simultaneous water shift and B 1 mapping (WASABI) approach of Schuenke et al 13 (t p = 5000 µs, B 1 = 3.7 µT, recovery time (saturated/M 0 ): 3 s/12 s, 31 equally spaced presaturation offsets from −1.8 to +1.8 ppm). This was executed with the same imaging parameters as the CEST MRI sequences. To enable correction for spatiotemporal changes in B 0 , 14 interleaved WASABI acquisitions were performed before and after each CEST examination.
For quantification of tSNR, 16 repeated measurements were performed without any presaturation module and a pause of 3 seconds in between consecutive image acquisitions. Different voxel sizes, acceleration factors, and nominal excitation FAs with otherwise identical parameters, unless stated differently, were investigated in preparation of the final 3D-EPI CEST protocol.

| CEST data evaluation
Data processing was performed using in-house written MatLab R2018a (The MathWorks, Inc., Natick, MA, USA) code, if not stated differently.
For both imaging methods, 3D-EPI and 3D-GRE, the postprocessing was identical. All images of the same examination were registered to the unsaturated image of the first WASABI to correct for subject movement during the acquisition. This was done automatically using the AFNI toolbox of Cox et al. 15 Afterward, the brain data were segmented automatically using SPM12 16 in case of the EPI and manually in case of the GRE readout to extract only the brain volume (including CSF, gray [GM], and white matter [WM]). The segmentation was of importance for the subsequently performed principal component analysis (PCA) based denoising. The CEST data sets were then corrected for dynamic changes in B 0 14 using ΔB 0 (r,t) maps derived from the interleaved WASABI measurements. To reduce noise in the images, PCA denoising 17 was applied next. Median criterion with correction factor β = 1.29 18 was used to determine the cutoff eigenvalue in the reconstruction. Malinowski indicator function was exemplarily applied to one data set as well. With the so-called real error RE, 19 the cutoff eigenvalue k ind can be derived as the minimum of RE(k)/(c-k) 2 , where c is the number of columns (rows are observations in Casorati notation; columns are reshaped 3D volumes at a single offset frequency) in the considered matrix. 20 The denoised data were then corrected for spatial B 1 inhomogeneity, including three different nominal B 1 values for Z-B 1 correction. 3 To quantify the CEST effects, a 4-pool Lorentzian model 21 was fitted to the postprocessed data voxel by voxel using a nonlinear least squares algorithm. The peak positions of APT, rNOE, and semisolid magentization transfer (ssMT) were fixed with respect to the frequency offset of direct water saturation (Supporting Information Table S1). To evaluate CEST effects within a certain tissue type, the tissue probability maps as derived with SPM12 were cut at a certain threshold (50% for GM and WM; 25% in case of CSF). Lower threshold for CSF was chosen to ensure contributions by CSF were separated from GM and WM.
For estimations on the point spread functions (PSF) of different readouts, numeric simulations had been performed for which a description can be found in the Supporting Information. These simulations were also used to investigate the effect of different FAs.
To evaluate the reproducibility of the derived results the Lorentzian amplitudes of APT, rNOE and ssMT pools were compared for different acquisitions. Within the tissue masks of GM and WM, the median amplitudes were determined. These were in the next step compared to the results of the other acquisitions in terms of mean value and SD, combined as a coefficient of variation (COV), which is the ratio of SD over mean value.

| tSNR evaluation
tSNR was determined in each voxel as the ratio of mean value over SD for repeated measurements. As shown in Figure 1A, the tSNR decreased on average by approximately 30% going from (2.0 mm) 3 to (1.8 mm) 3 isotropic resolution. The increased average tSNR for higher FAs went along with a less homogeneous spatial distribution ( Figure 1B,C and Supporting Information Figure S1). For FA > 20° a significant decrease of tSNR in the center and in lower regions of the brain was observed. For the final protocol, a FA of 15° was chosen. The PSF simulations ( Figure 2) revealed a relative FWHM 2 of on average 1.33 for FA = 15°. The maximum relative magnitude was found to be 4.4% for FA = 17.5° (Supporting information Figure S4A). These numbers also prove a FA of 15° to be suitable. Increasing the total acceleration factor from 4 to 6 decreased tSNR by around 17% (Supporting Information Figure S2), which approximately corresponds to the expected thermal noise increase, that is, the expected additional g-factor penalty seems to be compensated for by reduced physiological noise (ratio of square root of undersampling factors: 18%). The non-blipped sixfold CAIPIRINHA acceleration nearly preserved signal with similar echo train length (24.3 vs. 23.9 ms, same EPI factor of 32) and allowed fewer excitations for the complete readout (45 vs. 66). The latter likely reduced physiological contributions to the tSNR. 22 For the final protocol (FA = 15°, non-blipped CAIPIRINHA = 1 × 6 (shift = 2) , partial Fourier = 6/8, 1.8 mm isotropic resolution), the tSNR of the EPI readout was ≈ 75 even in the cerebellum. In contrast to the previously published protocol for 7T 4 , we worked out and optimized a modified 3D-EPI readout, for example, in terms of undersampling (GRAPPA vs. CAIPIRINHA), full versus elliptical sampling, and segmentation of one versus three. This yielded ≈20% higher tSNR for the novel 3T protocol (Supporting Information Figure S3). The PSF analysis revealed that the referring FWHM increased by 8.7% in the first and 9.5% in the second phase-encoding direction (compared to the 7T protocol) with maximum 1.75 pixels in the first PE direction ( Figure 2E,F). Compared to the unfiltered PSF, the two protocols yielded a relative FWHM of 132 and 143%, respectively, in the first phase-encoding direction and 112 and 123% in the second.

1
As shown in Figure 3A to C, the B 1 amplitude varied by approximately 50% across the whole brain volume. The deviation along the head-foot direction was larger than the deviation along the other two (left-right and anteroposterior) directions ( Figure 3A).
If the data were not corrected for B 1 inhomogeneity, fitted CEST contrasts were significantly altered as shown for rNOE in Figure 3D,E. For comparison of B 1 correction methods, the FWHM of Z-value distribution in GM and WM was determined by fitting a Gaussian distribution at each presaturation offset. In Supporting Information Figure S5A, it is shown that using as little as n = 2 B 1 values for correction already yielded 40/48% narrower FWHM in Z-value distributions in GM and WM than without B 1 correction. The resolution. Higher FAs provide higher average tSNR but at the cost of increased spatial heterogeneity (see also (C)), with very low tSNR in the center and lower regions of the brain. Additionally, a single image of the time series is shown for both FAs. (C) Closer look at FAs around 15 degrees with mean and SD in different axial slices. (D) Comparison to the GRE readout at 2.34 mm. EPI at 1.8 mm isotropic resolution outperforms the 2.34 mm GRE in terms of tSNR. In the cerebellum, the tSNR is around 80. FA, flip angle; GRE, gradient echo; tSNR, temporal SNR assumes linear behavior of Z-spectra over a 50% change in B 1 . If only n = 2 B 1 values were included, the fitted amplitudes were altered (Supporting Information Figure S5B). Supporting information Figure S5B to E also show that including an increasing number of B 1 values in the correction reduced noise in the fitted contrast maps. This held true for spline but not necessarily for linear interpolation. We therefore suggest to measure at three different B 1 amplitudes and to use spline interpolated Z-B 1 correction.
As shown in Figure 6, fitted CEST pool amplitudes revealed reproducible GM/WM contrast for different healthy subjects. Reproducibility was investigated both across the same subject at different sessions and for direct repetition within the same with CAIPIRINHA sampling in (C); and in (B,D) the resulting relative signal intensity is shown. PSFs (E,F) were determined from Fourier transformation of a 2D delta function that was convoluted with the signal evolution (B,D) determined from k space sampling taking into account both T 1 and T 2 (nominal matrix size (RO × PE × 3D): 144 × 126 × 88, FA = 15°, T 1 = 1300 ms, T * 2 = 45 ms and initial magnetization M 0 = 1 for both protocols. Variable TR for different excitations in 3T protocol with semielliptical sampling). Gamma is the FWHM of the PSF in both subplots. For the 3T protocol, like regular centric-out sampling, each 3D phase encoding slice starts with a new excitation; however, the first phase encoding is offset by ±1 line according to the CAIPIRINHA pattern. CAIPIRINHA, controlled aliasing in parallel imaging results in higher acceleration; PE, phase-encoding directions; PSF, point spread function; RO, readout; T, tesla session (intrasubject) as well as across different subjects (intersubject). The intrasubject coefficients of variation (COV) in fitted amplitudes were below 8.5% for examinations in the same session and below 7% for examinations in different sessions for APT, ssMT, and rNOE both in GM and WM ( Figure 7A,B). It was found that ssMT showed four times higher COV than the other fitted amplitudes when evaluated for examinations within same session. The intersubject COV was below 4% for n = 3 subjects ( Figure 6C) for all tissues and CEST pool amplitudes. In this case, none of the fitted amplitudes showed a manifold increased COV as compared to the others within both tissue types at the same time.

| Comparison with 3D-GRE protocol
Both GRE and EPI readout were compared at an isotropic resolution of 2.34 mm with same CEST presaturation module. Still, even at a higher resolution of 1.8 mm, the EPI (black dash-dotted line, Figure 1A) outperformed the GRE at 2.34 mm isotropic resolution (red dashed line, Figure 1D) in terms of tSNR. At (2.34 mm) 3 , the tSNR of EPI was approximately two times higher than that of the GRE ( Figure 1D). CEST measurements were acquired in the same volunteer but in different sessions. EPI was applied with matrix size 110 × 96 × 72 (RO × PE × 3D: head-foot × anteroposterior × left-right) and had a total acquisition time of 4:09 minutes for 58 offsets including the unsaturated M 0 image. On the contrary, the GRE was executed with matrix size 96 × 78 × 72 (RO × PE × 3D: head-foot × anteroposterior × left-right) and took 4:59 minutes. Readout durations were 1.7 seconds for GRE and 0.8 seconds for EPI, which also showed less blurring in the raw images (Figure 7, last row).
Considering the fitted CEST amplitude maps (Figure 7), it was found that the EPI readout allowed to detect finer anatomical structures in the brain. Especially the fitted amplitudes of rNOE were more blurry using the GRE readout. Still, both readouts yielded fitted CEST amplitudes in GM and WM segments that deviated by less than 10% ( Table 1). The contrast ratios between GM/WM were found to be comparable with {0.9, 1.2, 1.6} (EPI) and {0.9, 1.1, 1.4} (GRE) for APT, rNOE and ssMT.

| Influence of denoising
To estimate systematic effects of PCA denoising on the appearance of the Z-spectrum, the average spectrum in a homogeneous WM region of interest was considered. Supporting information Figure S6 shows that the residual Z-values due to PCA denoising according to Median criterion 17 −5.6 × 10 −6 , and the average spatial SD in the non-denoised Z-spectrum was 1.0 × 10 −3 . The resulting fitted amplitudes did not reveal anatomical structures depending on the denoising, as shown in Supporting Information Figure S7 and S8. In addition to the suggested postprocessing that includes PCA denoising using the Median criterion, stronger denoising with application of Malinowski's indicator function 20 was investigated. This had previously been applied in the studies of Goerke et al 24 and Deshmane et al. 6 In our study, Malinowski's indicator function suggested to include roughly 40% of the number of components suggested by the Median criterion. It was observed that this made, for example, the structures in the striatum-especially the putamen-more clearly visible in the fitted amplitude maps (Figure 8). The differences in fitted amplitudes between the two cutoff criteria did not show obvious anatomical structure besides some areas in the CSF (Supporting Information Figure S9B,C,D). The coefficients of variation of the medians of fitted amplitudes in GM and WM for both cutoff criteria were less than 1.4% for APT, rNOE, and ssMT (Supporting Information Figure S9A).

| DISCUSSION
It was shown that with the worked-out 3D-EPI readout, a tSNR of 75 could be achieved in all regions of the brain, , the approach at B 0 = 3T benefits from more homogeneous saturation and excitation, especially in lower regions of the brain. In addition, the 3D-EPI readout benefits from lower field strength such that T 2 is increased by about 65% compared to 7T. 11 Assuming T * 2 scales similar as T 2 with field strength, at the given echo train length of ≈ 25 ms this is 30% higher residual transversal magnetization at the end of the readout (assuming monoexponential decay and linear scaling of magnetization with field strength), compensating partly for lower SNR as compared to ultrahigh field. The final protocol based on 3D-EPI allowed to perform CEST MRI at 1.8 mm isotropic resolution within 4.3 s per presaturation offset (1.2 s readout duration) with 256 × 224 × 156 mm 3 whole brain coverage. A previous approach of Zhu et al 25 27 had chosen a FLASH approach with two shots per readout volume at B 0 = 7T. This enabled high in-plane resolution of (0.6 mm) 2 with a FOV of 140 × 140 × 48 mm 3 for 16 slices but took 16 s per saturation offset and did not provide whole brain coverage. Deshmane et al 6 worked with a 3D-GRE based snapshot approach, yielding (1.7 mm) 2 in-plane resolution for 18 slices with FOV = 220 × 180 × 54 mm 3 and 2.9 s readout. In general, snapshotbased approaches benefit from additional freedom in the presaturation such that no specific magnetization preparation has to be achieved for the image readout. It can be seen that our suggested protocol offers one of the highest spatial resolutions combined with a large FOV, it is significantly faster than most of the reported approaches and is flexible due to its snapshot realization.
One remaining drawback that all whole brain approaches share is the increased B 1 inhomogeneity within the large FOV even at clinical field strength. It has been shown previously 3,28 that an effective B 1 correction strategy is a necessary prerequisite for reliable CEST quantification at ultra-high field. However, if homogeneity over the whole brain volume is aimed at, then B 1 correction becomes necessary even at 3T, as shown in this study. This prolongs the scan time significantly because it demands repeated CEST acquisitions with different B 1 scaling. Similar to the ultra-high field study, 4 in the present study measurements at [0.75, 1.00, 1.25] times the nominal B 1 were shown to

F I G U R E 6 Comparison of fitted CEST contrasts in GM and WM of healthy subjects. Within same session and volunteer (A) and in different sessions (B). Across n = 3 sessions and volunteers (C)
; variations in fitted CEST contrasts calculated as (SD/mean over tissue type) are comparable within the same volunteer and across volunteers. Fitted CEST contrasts are stable with coefficient of variation (SD/mean) less than 9%. APT, amide proton transfer; GM, gray matter; rNOE, relayed nuclear Overhauser enhancement; ssMT, semisolid magnetization transfer; WM, white matter be sufficient for spline interpolated Z-B 1 correction. The range of nominal B 1 values covered 95% of all occurring B 1 values and avoided extrapolation of the Z-spectra when correcting for B 1 inhomogeneity. The FWHM of the distributions of Z-values for all offsets were investigated in GM and WM separately after Z-B 1 correction. Including the three different B 1 values mentioned above and performing spline interpolation gave smaller FWHM than linear interpolation from five different B 1 values. On the contrary, including up to five B 1 values for spline interpolation only marginally decreased the FWHM further, indicating that the width of the natural distribution within the region of interest was already reached. Correction using two B 1 values yielded 44% smaller FWHM of Z-values than without correction. Still, in this case the fitted Lorentzian amplitudes were visibly altered. If parallel transmit systems are available, B 1 inhomogeneity might be further mitigated using B 1 shimming 29 or the MIMOSA approach, 30 which could directly be combined with our protocol.
If only certain brain areas are of interest, the number of scans at different presaturation B 1 values may be reduced given by the B 1 inhomogeneity in the considered slices. For example, for a single, more cranially located axial slice in which also B 0 distortions due to the air cavity of the mouth are less dominant, the B 1 deviation was below 20% ( Figure 3A, Supporting Information Figure S10).
We also suggest to perform interleaved B 0 measurements after each CEST scan. This might not be of major importance as long as peaks can be fitted to the Z-spectrum. But in case that high power presaturation modules are used or asymmetry analysis is performed for evaluation, 31 alteration of B 0 inhomogeneity over time should be taken into account. 14 Necessity of repeated acquisitions at different presaturation B 1 also increased examination time, and therefore, the sensitivity to temporal variations in B 0 inhomogeneity. With accelerated field mapping methods, for example, DREAM/3DREAM 32,33 for B 1 and dual-echo EPI for B 0 mapping instead of WASABI, the duration of field mapping F I G U R E 7 Fitted amplitudes of CEST pools: APT, rNOE, and ssMT. 2.34 mm isotropic resolution GRE (left) and EPI (right) both with same CEST presaturation module, same postprocessing, and in the same volunteer. Different segments were chosen for GRE and EPI readout, respectively, because data were not coregistered onto each other. Both readouts reveal similar GM and WM contrast in the fitted amplitudes of APT, rNOE, and ssMT. With the EPI readout, more detailed anatomical structures are visible. Additionally, in the last row unsaturated images are shown for both readouts. APT, amide proton transfer; GM, gray matter; GRE, gradient echo; rNOE, relayed nuclear Overhauser enhancement; ssMT, semisolid magnetization transfer; WM, white matter could be strongly reduced. Assuming this enables field mapping within one minute, the total duration for the suggested procedure would be reduced by 25%. Although shorter scan durations may result in less intervolume head motion over time, retrospective motion correction before CEST analysis is still required (see Supporting Information Figure S12). Strong imaging acceleration is also an effective means to reduce intravolume motion sensitivity. Here, whole-brain k-space acquisition in only 1.2 s over 45 excitations resulted in no apparent motion artifacts. However, this may have to be reevaluated for motion-prone subject groups.
To quantify CEST effects, the amplitudes of the fitted 4-pool Lorentzian model 21 were considered. Comparable fit models had recently been applied by, for example, Deshmane et al, 6 Goerke et al, 24 and Akbey et al. 4 Medians of fitted amplitudes in GM and WM were determined separately for each subject (Supporting Information Table S2). In Table 1, the mean values over different subjects are listed. In Figures 4, 5, and 8, as well as Supporting Information Figure S13, it can be observed that especially the APT maps show some areas of strongly increased amplitudes in all healthy subjects. These match the expected spatial distribution of vessels known to have different CEST properties compared to brain tissue, 34,35 which might explain the observations in APT maps. More data on this issue can be found in Supporting Information Figure S13 to S15. Comparing the fitted CEST amplitudes to the results of Deshmane et al 6 (B 0 = 3T, 3D-GRE readout) showed similar behavior in GM relative to WM for APT, rNOE, and ssMT (Table 1). Nonetheless, absolute values of fitted amplitudes differed. For APT and rNOE, we derived about twice as high amplitudes as reported by Deshmane et al. 6 For ssMT, we found 1.4 times smaller amplitudes. This can be explained by the use of inversion pulses in the other study, which are known to generate strong MT effects, 36 whereas we used longer pulses for the sake of spectral selectivity. Additionally, a different number of pools was fitted to the data in both studies. Recently, Akbey et al 4 presented a CEST study at ultra-high field (B 0 = 7T, 3D-EPI readout). The reported GM/WM ratios (MTR LD metric 23 ; F I G U R E 8 Effect of denoising on the subsequently fitted CEST pool amplitudes. APT, rNOE, and ssMT amplitudes fitted after postprocessing that included denoising according to Median (left) and Malinowski (right) criterion. Arrows highlight areas were denoising causes substantial differences in the spatial distribution of fitted amplitudes. Especially in the striatum, structures that might be of similar shape as caudate nucleus or putamen become visible when denoising with the more aggressive Malinowski criterion. APT, amide proton transfer; rNOE, relayed nuclear Overhauser enhancement; ssMT, semisolid magnetization transfer data extracted from Figure 9 in Ref. 4) of APT, rNOE, and ssMT {1.1, 0.9, 0.8} were similar to what we observed at B 0 = 3T. Nevertheless, it becomes obvious at this point that for all different presaturation and Lorentzian models used, standardization is missing. Most importantly, the compared EPI and GRE readouts yielded the same contrast for the same CEST presaturation and postprocessing used (Supporting Information Figure S14).
To investigate reproducibility of the suggested CEST pipeline, the median of the fitted amplitudes was considered. The COV across three subjects and sessions was below 4% for APT, rNOE, and ssMT. Still, it was up to 8.3% for ssMT in GM when comparing the results of repeated examinations of the same subject in the same session. On the contrary, for APT and rNOE amplitudes, it was below 2% for these two data sets. The increased deviation in the ssMT amplitudes may have several contributions. One explanation could be different GM and WM segmentation that was performed independently for different examinations. Segmentation tissue volume of both data sets differed by 0.45% for WM and 0.87% for GM. Because ssMT amplitudes show larger GM/WM contrast than those of APT and rNOE, a different segmentation should affect ssMT stronger. Still, segmentation could not exclusively be concluded to explain the observed deviations in the ssMT amplitudes. The deviations rather indicate that with the suggested CEST pipeline fitted amplitudes can at worst be reproduced with a COV of less than 8.5%.
The determined Lorentzian peaks might be used to calculate more sophisticated CEST contrasts such as the apparent exchange-dependent relaxation, 37 which additionally considers differences in T 1 relaxation times. One exemplary apparent exchange-dependent relaxation evaluation is shown in the Supporting Information Figure S11, with T 1 estimated using the proposed 3D-EPI readout in a saturation recovery protocol.
We also directly compared the proposed 3D-EPI protocol to an established 3D-GRE-based whole brain protocol. 12 EPI outperformed the GRE readout in terms of tSNR by a factor of two at the highest possible nominal resolution of the GRE, which was 2.34 mm isotropic. In addition, for the described settings the EPI readout was more than two times faster than the GRE. This gain in tSNR, along with the reduction in measurement time, allowed increased spatiotemporal resolution in the CEST acquisition as presented in this work. Fitted CEST amplitudes revealed much more detailed anatomical structures of the brain in case of the EPI readout, as shown in Figure 7. Still, GM/WM ratios of fitted amplitudes derived with the established GRE readout {1.06, 0.89, 0.71} were comparable to that observed with EPI (Table 1). It seems that actually the 3D-GRE was at its limits in case of the presented whole brain settings, for example, leading to more blurry images compared to the 3D-EPI and an established slab-selective 3D-GRE. The latter revealed contrast maps very similar to the nonselective 3D-EPI (Supporting Information Figure S14).
CEST spectra can be strongly affected by direct saturation of fat signals. 38,39 Thus, a water excitation enables mitigation of this artifact. However, binominal water excitation pulses 9 are rather long (additional ≈ 1.1 ms per excitation for both GRE and EPI) and would prolong a GRE readout significantly. Fortunately, for the EPI much fewer excitations are needed compared to GRE (45 at 1.8 mm isotropic resolution vs. 487 for the GRE at 2.34 mm), and these longer pulses can easily be afforded regarding scan time for EPI. Thus, artifacts in the Z-spectrum originating from fat signals are directly suppressed in snapshot EPI. In general, the faster readout of EPI is beneficial when comparing the PSF of GRE and EPI. Overall T 1 signal decay during readout can be described by a Look-Locker decay 2,40,41 in both cases. Assuming a T 1 of 950 ms, 11 this yields a decay to 6.4% of the initial signal amplitude for EPI after acquisition of one 3D volume and 0.55% for GRE (EPI at 1.8 mm and GRE at 2.34 mm isotropic resolution). For EPI along the first PE direction (anteroposterior), signal decay is governed by T * 2 . Still, with T * 2 = 30 ms, approximately 40% of the initial signal remains at the end of the considered PE-line (one PE-line containing up to 32 readout lines was encoded per excitation).
To improve SNR, PCA denoising 17 was applied during the postprocessing for both GRE and EPI, with cutoff eigenvalue determined according to Median criterion. A more aggressive denoising as introduced by Malinowski 20 was found to provide better visibility of anatomical structures in the striatum. We found that the COV due to different denoising was below 1.4% for the fitted Lorentzian amplitudes. Difference maps of amplitudes revealed no anatomical structures, but contiguous areas with larger deviations that are not noise-like (Supporting Information Figure S9) were observed. Also, Breitling et al 17 reported that Malinowski indicator function can be problematic if applied to determine the cutoff eigenvalue for CEST data. In accordance with their findings, we therefore suggest to apply the more conservative Median criterion. The chosen denoising yielded noise-like residuals when considering a region of interest in WM with and without denosing and did not alter the fitted amplitudes (Supporting Information Figure S6 to S8), which provides strong evidence for the assumption that the denoising did indeed only remove noise but did neither remove nor add systematic features to the Z-spectra.

| CONCLUSION
It was demonstrated that it is possible to acquire reproducible (coefficient of variation of less than 9%) CEST contrasts with whole brain coverage at B 0 = 3T in 4:22 minutes per presaturation B 1 . In addition, with 1.8 mm isotropic resolution we were able to decrease the voxel size by 60% compared to currently used multislice or 3D CEST protocols. The suggested protocol is adaptable and may be used for any CEST presaturation module and fit model of interest. It therefore increases the opportunities of CEST investigations at B 0 = 3T widely.

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the Supporting Information section. Influence of denoising according to the Median criterion on the mean Z-spectrum in a white matter region of interest in a healthy volunteer. The mean Z-spectrum did not significantly differ after denoising which is shown, for example, in the difference of the Z-spectra. The mean value of difference is −5.6 × 10 −6 , whereas the spatial standard deviations are on average 1.0 × 10 −3max (without PCA) and 9.4 × 10 −4 (including PCA), respectively FIGURE S7 Comparison of fit results in a healthy volunteer with respect to the effect of PCA denoising (cut off eigenvalue determined according to the Median criterion). Plots show correlation of fitted Lorentzian amplitudes with and without PCA denoising. Legends provide: Pearson r, P-value as well as slope and offset of a linear fit to the data. High correlation of r > 0.92 for all contrasts (APT, rNOE, ssMT and (direct saturation) DS) indicates that PCA denoising does not intorduce any systematical bias to the fit results FIGURE S8 Fitted Lorentzian amplitudes of APT, rNOE, SSMT and DS (direct saturation) are shown in the same axial slice separately with and without PCA denoising (cut off eigenvalue determined according to the Median criterion). Last row shows the difference between the two sets of fitted amplitudes. The difference maps do not reveal any significant anatomical structure that would have potentially been introduced due to the PCA denosing region of the brain that had more homogeneous B 1 distribution than other areas of the brain. B 1 correction including only two different B 1 values (linear interpolation) might, therefore, be sufficient. Still including five B 1 values and applying spline interpolation yields smoother maps due to averaging FIGURE S11 CEST MRI acquired according to the final protocol combined with a saturation recovery sequence to determine T 1 and from this calculate the T 1 corrected AREX (37) contrast FIGURE S12 Same data and post processing but once with and once without motion correction. Arrows highlight areas were substantial differences in the fitted Lorentzian amplitudes can be seen if motion correction is not applied. Even though the optimized protocol allows relatively fast CEST MRI acquisition motion would significantly corrupt the results if not corrected FIGURE S13 Fitted Lorentzian amplitudes (rNOE, APT, ssMT and (DS) direct saturation) of four healthy volunteers. For all volunteers only in the APT maps voxels with much larger amplitudes are visible FIGURE S14 Fitted Lorentzian amplitudes and maximum intensity projections for an established slab selective 3D-GRE (2,6) and the novel non-selective 3D-EPI readout. In both cases only in the APT maps similarly shaped and located areas show higher amplitudes (see arrows). Both readouts were acquired in different volunteers, and therefore, have different slice orientations. Data includes slices 2 to 11 (GRE) and 87 to 100 (EPI), respectively FIGURE S15 Same 3D-EPI measurement as shown in Supporting Information Figure S15. Residual sum of saqures (RSS) was determined in slice 87 to 100. It is found that RSS does not increase with increasing APT amplitudes (A) and that highest RSS is observed in the ventricles TABLE S1 Fit parameters for all i = 0, …, 3 amplitudes (Ai), peak widths (Gi), initial magnetization (Zi) and position of the bulk water pool (dwi). Positions for i = 1, …, 3 were fixed at 3.1, −1.25 and −3.1 ppm relative to dwi. Other parameters fitted with lower bond (lb), upper bond (ub) and start values (p0) TABLE S2 Resulting amplitudes of 4-pool Lorentz fit for n = 5 different examinations in 3 different subjects (S0, S1 and S2). All values are given in percent to increase visibility. For each examination and subject median and standard deviation (std) were determined in gray (GM) and white matter (WM) for APT, rNOE and ssMT. In