Enhanced and Robust Contrast in CEST MRI: Saturation Pulse Shape Design via Optimal Control

Purpose: To employ optimal control for the numerical design of CEST saturation pulses to maximize contrast and stability against $B_0$ inhomogeneities. Theory and Methods: We applied an optimal control framework for the design pulse shapes for CEST saturation pulse trains. The cost functional minimized both the pulse energy and the discrepancy between the corresponding CEST spectrum and the target spectrum based on a continuous RF pulse. The optimization is subject to hardware limitations. In measurements on a 7 T preclinical scanner, the optimal control pulses were compared to continuous-wave and Gaussian saturation methods. We conducted a comparison of the optimal control pulses were compared to with Gaussian, block pulse trains, and adiabatic spin-lock pulses. Results: The optimal control pulse train demonstrated saturation levels comparable to continuous-wave saturation and surpassed Gaussian saturation by up to 50 \% in phantom measurements. In phantom measurements at 3 T the optimized pulses not only showcased the highest CEST contrast, but also the highest stability against field inhomogeneities. In contrast, block pulse saturation resulted in severe artifacts. Dynamic Bloch-McConnell simulations were employed to identify the source of these artifacts, and underscore the $B_0$ robustness of the optimized pulses. Conclusion: In this work, it was shown that a substantial improvement in pulsed saturation CEST imaging can be achieved by using Optimal Control design principles. It is possible to overcome the sensitivity of saturation to B0 inhomogeneities while achieving CEST contrast close to continuous wave saturation.


Introduction
Chemical Exchange Saturation Transfer (CEST) is an emerging MRI technique that allows for detection of molecules via protons that are in chemical exchange with the surrounding water protons.This exchange transfers saturation from molecules to water, which can be detected indirectly by standard MRI methods.This enables fast and high-resolution images that contain spectroscopic information in every voxel.In CEST imaging, the signal obtained from a mobile proton pool gets amplified, which allows for the detection of low concentration molecules [1,2].Contrast of relevance for clinical applications can be derived from groups such as amide, amine, and hydroxyl in mobile protons and peptides, as well as from glycosaminoglycans, glutamate, and shift in creatine and pH levels [3][4][5][6][7][8][9].Early uses of CEST MRI in medical environments have shown promise for evaluating various medical conditions, encompassing strokes, cancer, osteoarthritis, muscle function, lymphedema, multiple sclerosis, Alzheimer's disease and other neurological pathologies.[10,11] A typical CEST experiment consists of the application of an off-resonant saturation pulses with a duration of a few seconds, followed by an imaging sequence.This is repeated for various off-resonant frequencies to accumulate a CEST-or Z-spectrum.The quality of the Z-spectrum is not only crucial for the acquisition of accurate spectroscopic information, but also for the reliability and resilience of the resulting image data.One of the most important parameters influencing the characteristics of the Z-spectrum is the saturation strategy.A common approach in preclinical research is continuous wave (CW) saturation, which result in a spectrum with high exchange weighting.However, on clinical systems, CW pulses with a duration of several seconds cannot be applied due to hardware limitations.Consequently, pulsed CEST saturation is the standard on clinical scanners.[12][13][14][15].
The optimization of RF pulse parameters within a pulsed CEST experiment to increase CEST contrast was explored in the past [13,16,17].Nevertheless, a large number of pulsed CEST applications still rely on predefined saturation pulse shapes, such as Gaussian or a slight modification of Gaussian line shapes as well as Block Pulse Trains (BPT) [18].The research conducted by Yoshimaru et al. indicates that the shape of the saturation pulse is a critical parameter to consider when striving for optimal CEST spectrum quality [19].With the off-resonance Spin Lock (SL) method, an improved pulsed saturation technique was introduced [20].To reduce the sensitivities to field inhomogeneities, the tip-down and tip-up pulses were substituted with optimized adiabatic pulses for certain B 0 field strengths [21,22].RF pulse design by Optimal Control (OC) has been used in the past to design RF pulses within a constrained parameter space.By accounting for physiological and technical restrictions, this approach ensures the feasibility of RF pulses on clinical scanners, while minimizing certain parameters such as pulse duration, specific absorption rate or B 1 /B 0 sensitivity [23][24][25][26][27][28].
In this work, we design the RF pulse shape of a CEST saturation pulse train by OC.The design criteria were chosen to attain high quality saturation, including high exchange weighting and stability against field inhomogeneities.We compare the OC pulses with CW saturation on a 7T preclinical scanner in phantom measurements.We further compare our pulses with Gaussian, BPT and adiabatic SL pulses on 3T in phantom measurements.Moreover, we perform Bloch-McConnell simulations of the magnetization vector dynamics with OC and BPT to understand their behavior under the influence of B 0 inhomogeneities.

RF pulse optimization and Bloch-McConnell simulation
The RF magnitude, B 1 (t), of a pulse train with a saturation time of T sat , duration of a single pulse t d and pause between pulses t p is optimized by minimizing the objective function, min The first term of the cost functional minimizes the pulse energy with regularization parameter α > 0. While M z (ω, T sat ), M zdes (ω) denote the spectra at the end of the pulse with two pools, Mz (ω, T sat ), Mzdes (ω) represent the spectra with solely water pool.A is the Bloch-McConnell system matrix depending on RF pulse B 1 (t).M (ω, t) and M (ω, t) are the solutions of the Bloch- McConnell equations for the off-resonance ω in the z-spectrum Ω.The simulation is discretized in time t with a step size of ∆t = 100 µs.Temporal discretization was chosen based on [30].The second and third terms minimize the difference between the CEST spectrum, M z , at the end of the RF pulse and a desired target CW CEST spectrum, M zdes .σ 1,2 > 0 are weighting parameters.The metric for the differences is conducted via the L p -norm, p ≥ 2 even.ǫ specifies the allowed deviation from the target magnetization.We loop over p starting with p = 2 increasing the value successively as proposed by [27], since a lower value of p has been found to confer advantages for globalization.
Optimization was constrained to the points during t d .The pause t p was set to B 1 (t) = 0.
The numerical optimization itself is based on a trust-region semi-smooth quasi-Newton method [27].
Exact discrete derivatives are supplied by adjoint calculus [29].First, the adjoint equations to the Bloch-McConnell equations are introduced as ). [3] Here, q is the adjoint state corresponding to M .Analogously the other adjoint state q can be derived with M , in which only the water pool is active.Subsequently, the reduced gradient g of the cost functional J can be identified as whereby The CEST-spectra are simulated using two-pool Bloch-McConnell equations solved numerically by symmetric operator splitting [30].
For simulation and optimization, the water pool was set to relaxation times T 1w = 1000 ms and T 2w = 80 ms.The solute pool was simulated to resemble a creatine phantom with relaxation times T 1s = 1000 ms and T 2s = 160 ms and a fraction rate of f s = 0.002 [31].The exchange between the solute pool and the water pool was set to k sw = 210 s −1 .The exchange was estimated for creatine at room temperature and pH 7.4 with the formula suggested in [33].The B 0 field was 3 T.The spectrum was optimized for Ω ∈ [−6, 6] ppm, discretized in ∆ω = 0.02 ppm steps with a total number of spatial points of N ω = 601.M zdes was simulated with the same parameters for a T sat = 1 s CW saturation pulse with an RF amplitude of B 1 = 1 µT.The optimized pulse was constrained to a saturation time of T sat = 1 s and had a duty cycle of 90 %, a pulse duration time of t d = 100 ms, and t p = 12.5 ms pauses between each of the nine pulses (Pulses optimized with different DC, T sat , t d and t p can be found in the Supporting Information Figures S5-8).Optimization was constrained to the points during t d .The pause t p was set to B 1 (t) = 0.The RF amplitude was constrained to 0 ≤ B 1 ≤ 3 µT.The regularization was α = 100 σ 1 and σ 1 = σ 2 .We loop over p starting with p = 2 increasing the value successively as proposed by [27], since a lower value of p has been found to confer advantages for globalization.The optimization was initialized with a BPT and converged with p = 4.

Evaluation of OC pulse in simulations
We evaluated the impact of various parameters on the saturation of the OC pulse in simulations, such as T 1 , T 2 , exchange rate k, B 0 and B 1 amplitude scaling of the pulse.The performance of the OC pulse was assessed by the smoothness of the spectra, the stability of the Magnetization Transfer Ratio Asymmetry, M T R asym , and the maximum value of M T R asym across a range of parameters.
To assess the impact of Magnetization Transfer (MT) on the performance of the OC pulses in comparison to Gaussian and BPT saturation, simulations were carried out using a three-pool model at 3 T.The MT pool was based on the WM parameters described in [32].

Experimental setup
Performance of OC saturation was compared to a CW pulse and Gaussian pulse train.These pulses were implemented on a 7 T preclinical scanner system (Bruker BioSpec 70/20 USR, Ettlingen, Germany).The CW pulse was applied with T sat = 1 s and B Phantom measurements were performed on a Siemens Vida 3 T clinical scanner system (Siemens Healthineers, Erlangen, Germany) with a 20 channel head coil.OC, SL, Gauss and BPT were implemented in Pulseq-CEST (Pulseq release 1.3.1)[35,36].The saturation time was set to T sat = 1 s, with a Duty Cycle (DC) of 90 %.The B 1,RM S was set to 1 µT, this means that all pulse trains had the same energy and therefore the same specific absorption rate (SAR).The OC and Gauss train were implemented using a constant phase.The phase of the BPT was adjusted after each pause by ω n T pulse B 0 γ, with the current offset ω n , to ensure a smooth spectrum of the BPT.The adiabatic SL pulse was provided by the Pulseq-CEST library.This pulse was optimized for 3 T as described by Kai Herz et al. in [22].The adiabatic tilting pulses were factored into the calculation of the pulse block duration, as well as the B 1,RM S .The Gaussian saturation as it appears in Pulseq-CEST is windowed with a cosine.(For more details on the pulse trains used, see the Supporting Information Figure S1.) The pulse trains were applied to a 87 mM creatine monohydrate phantom in a 500-ml glass sphere.
The relaxation times were reduced with Manganese(II)chloride to T 1 = 680 ms, T 2 = 90 ms [37].The pH was stabilized at 7.0 with a phosphate buffer.Furthermore, an air bubble was placed intentionally at the top of the glass sphere to introduce B 0 inhomogeneities.CEST-spectra were measured for ω ∈ [−6, 6] ppm with a step size of dω = ∆ω = 0.1 ppm.For readout a centric reordered 2D-GRE was used, with parameters: F OV = 128 × 128 mm 2 , flip angle α = 10 • , T R = 12 ms, T E = 40 ms, base resolution res = 64×64 and slice thickness 5 mm with one slice through the center of the sphere.
In addition, to investigate the influences of the B + 1 field on the acquired Z-spectra, a B + 1 -map was estimated from the Bloch-Siegert (BS) shift [38] using a 3D-GRE readout.The sequence parameters were as follows: T R = 150 ms, α = 30 • , matrix size = 64 × 64 × 64, FOV = 128 × 128 × 128 mm 3 , BS pulse with an off-resonance frequency of 4 kHz and a nominal BS-shift of φ BS = 0.78 rad.The B + 1 -maps were denoised using a 5 px median filter.

Postprocessing and analysis of experimental data
To evaluate the exchange weighting of the different pulses, M T R asym [4,39] was used.A pixel-wise two pool Lorentzian fit using the MATLAB function lsqcurvef it was employed to assess and correct the B 0 shift and to determine the peak of the CEST effect [31,40,41].To improve the visualization of structural inhomogeneity, the CEST images were denoised using TGV.[42,43].
The geometry of the spherical glass phantom with an air bubble on top introduced different levels of B + 1 and B 0 inhomogeneities into the measurement.The frequency shift in the B 0 map was calculated in every pixel with the water peak shift of the fitted Gauss saturation.Spectra were evaluated in different ROIs placed over the phantom to investigate different combinations of B + 1 and B 0 .Moreover, the robustness of the spectra and the quality of the corresponding fit were examined by inspecting the Z-spectrum of representative pixels in the inhomogeneous region.

Simulation of magnetization vector dynamics of OC pulses
To investigate and visualize the dynamic behavior of the magnetization during saturation with OC and BPT, the x-, y-and z-components of the magnetization were simulated over the duration of the pulse train in the water pool at -1 ppm.The simulation includes two scenarios: one in which the phase remains constant throughout the duration of the pulse, as used during the optimization process, and another in which the phase was adjusted after each pulse by ω n T pulse B 0 γ.The latter scenario is necessary for obtaining a smooth BPT spectrum.Alterations in the phase of a pulse can be seen analogously to changes in either the pulse center frequency ω n or the B 0 field (see Supporting Information for details).The OC pulse yielded a smooth spectrum when T 2 ≤ 0.15 s.However, when T 2 was increased to 0.25 s, the spectrum began to exhibit wiggles.In (c) OC saturation was simulated for exchange rates k ∈ {50, 100, 200, 500, 1000, 2000} s −1 .The OC pulse produced a smooth spectrum for all values of k, however the M T R asym peak showed minor wiggles for k = 2000 s −1 .In (d), the OC pulse showed a smooth spectrum for all investigated field strengths B 0 ∈ {1.5, 3, 5, 7} T. In (e), the OC pulse was scaled by a factor so that the B 1RM S value was B 1RM S ∈ {0.3, 0.5, 0.7, 1, 1.2, 1.4, 2} µT.When the scaling surpassed 200 %, or when the B 1RM S fell below 30 %, the saturation began to introduce minor distortions in the spectrum.The peak value of the OC M T R asym demonstrated a performance comparable to that of CW saturation across all parameters examined.The maximum of M T R asym consistently exceeded that of Gaussian saturation at every simulated point.

Pulse optimization and simulation for different parameter ranges
Simulations using a three-pool model, inclusive of an MT pool, demonstrated a general decrease in M T R asym across CW, OC, BPT, and Gaussian saturation.Relative to the two-pool model, the MT pool led to reductions in M T R asym of 40.43 % for CW, 43.97 % for OC, 42.43 % for BPT, and 43.93 % for Gaussian saturation (see Supporting Information Figure S9).

Comparison 7 T preclinical phantom measurements with simulation
Figure 2 compares results of the simulations with measurements at a field strength of 7T, for the CW, OC and Gauss pulse train.The CW pulse measurement shows a maximum M T R asym of 9.6 % and the simulation result matches with a maximum of 9.5 %.Similarly, the OC pulse measurement exhibits a maximum M T R asym of 9.1 %, while the simulation result shows a maximum of 9.0 %.
The Gauss pulse train measurement demonstrates a maximum M T R asym of 6.7 %, whereas the simulation result shows a maximum of 6.1 %.The results of these measurements indicate that the exchange weighting of the CW pulse is 4.4 % higher than that of the OC pulse and 40.3 % higher than that of the Gauss pulse train.Analogously, the simulation results show that the exchange weighting of the CW pulse is 6.7 % higher than that of the OC pulse and 57.4 % higher than that of the Gauss pulse train.

Creatine phantom measurements on a 3 T clinical scanner
Figure 3 shows the results of the 3 T phantom measurements for the OC, BPT, adia SL and the Gaussian saturation pulses.(a) depicts the uncorrected CEST image at the creatine peak of the M T R asym .For the different saturation strategies, the B 0 inhomogeneities had different effects on the images.Especially the BPT exhibited a high CEST peak shift and expressed wavelike behavior.
The shift is particularly expressed under the air bubble and on the edge of the sphere.Across all images, the center of the sphere exhibited a relatively low frequency shift and high homogeneity.
The OC, adia SL, and the Gauss saturation pulses exhibited a comparably homogeneous CEST contrast and did not express wave-like behavior in the images.After B 0 correction (b) and TGV denoising (c), the BPT still showed significant wave-like behavior at the largest B 0 offsets.The the images generated using the OC, adia SL and Gauss saturation pulses did not express these behavior.The OC pulse generated a homogeneous CEST image with high contrast, while the adia SL pulse and Gaussian pulse still displayed a slight inhomogeneous pattern.Specifically, at the bottom of the sphere, the adia SL pulse appeared to have a reduced CEST amplitude, whereas the BPT exhibited again wave-like behavior.The homogeneity increased after fitting the OC and the BPT images (d, e).In the BPT image, the fit was able to partially reduced the wavelike artifacts compared to the M T R asym images.The fitted adia SL image exhibited a slight degradation in homogeneity.Additionally, the fit indicated a higher exchange weighting for the Gauss pulse train, particularly at the center of the sphere.The two-pool fit showed lower noise overall.
The difference in exchange weighting and homogeneity among the different saturation regimes, particularly in the region of the highest B 0 shift, is shown by a vertical profile through the center of the sphere of the B 0 corrected TGV denoised M T R asym images (f).The wave-like behavior of the BPT is evident between pixels 10 and 30, with a slight peak in exchange weighting for the OC and adia SL and a comparable smooth Gauss saturation in this region.The exchange weighting of the Gauss pulse train is substantially lower compared to other saturation strategies, with the OC pulse exhibiting the highest exchange weighting, followed by the BPT pulse and the adia SL pulse.

Influence of field inhomogeneities to the spectrum
Figure 4 (a) illustrates the changes in ∆ω due to ∆B 0 .In (b), the B + 1 map divided by the nominal BS B 1 level, B 1BS , is presented along with a line profile.The profile was selected to traverse from a region with least ∆ω shift to one with high ∆ω shift, while also passing through the sphere's center to span the entire B + 1 range.Therefore, in one half of the profile, the influence of B 1 is isolated, while in the other half, both B 1 and B 0 field imperfections are present.To better visualize the impact of field inhomogeneities on the saturation method, the M T R asym images were normalized and scaled using the maximum value of 1 and the ∆ω as well as the profile in the scaled M T R asym images can be seen in (c).The B + 1 achieved its peak in the center of the sphere, reaching 108 % of the nominal B 1 level, and decreased to the lowest point, 91 %, at the edge of the sphere.The ∆ω remained low, < 0.3 ppm, until the center, then increased to a peak of 0.12 ppm.Within the homogeneous B 0 region, the CEST contrast was comparable for OC, adia SL, and Gaussian saturation, In comparison, the BPT M T R asym more closely aligned with changes in B + 1 .As the ∆ω shift increased, notable deviations emerged in the BPT.Similarly, the adia SL displayed deviations relative to the profile's low ∆ω region.In contrast, the OC and Gaussian saturation remained relatively stable even amidst increasing frequency shifts.Notably, the OC exhibited the highest and most consistent saturation, expressing a saturation comparable to the low ∆B 0 region and demonstrating a uniform saturation over the profile.
The Regions of Interest (ROIs) studied, each exhibiting varying degrees of ∆ω shift, are depicted in Figure 4 (a).The corresponding statistics can be seen in the boxplots (d) and in Table 1.
Particularly pronounced is the ∆ω shift beneath the bubble at the sphere's apex, with additional inhomogeneities identified at the base and sides of the sphere, intersected by regions of minor shifts.
In the most homogeneous ROI, ROI 1, the frequency shift was measured as ∆ω = (0.037 ± 0.009) ppm.Here, the OC pulse displayed the highest exchange weighting with a peak M T R asym of 16.8 %, followed by the BPT at 16.1 %, adia SL at 15.4 % and Gaussian saturation at 11.3 %.The OC had the lowest Interquartile Range (IQR) at 0.6 %.Conversely, ROI 2 represented the most inhomogeneous region with a substantial frequency shift of ∆ω = (0.29 ± 0.04) ppm.Again, OC exhibited the highest mean M T R asym peak, 14.1 %, followed by the adia SL at 13.7 %.Notably, BPT (11 %) and Gaussian (9.9 %) peaks substantially declined compared to those in the center, while the IQR varied significantly between OC saturation (0.9 %) and BPT (4 %).ROI 3, located at the edge of the sphere, displayed an intermediate frequency shift (∆ω = (0.063 ± 0.034) ppm) and greater B 1 deviation compared to the center.The mean peak was again highest for OC saturation (15.9 %), with BPT (15.0 %) and adia SL (14.4 %) following closely, and was the lowest with Gaussian saturation (10.8 %).ROI 4, sharing a similar B 1 deviation to ROI 3 but presenting the lowest ∆B 0 of ∆ω = (-0.010± 0.011) ppm, again OC contrast was highest (16.0 %), with BPT (15.5 %), SL (15.1 %) and Gaussian (10.6 %) in succession.In this ROI the IQR was highest for the adia SL and comparable for the other saturation pulses.
Figure 4 (e) illustrates the mean spectra within ROI 2. The adia SL presents asymmetric saturation around 0 ppm and a plateau between 0 and 0.5 ppm, showcased in the uncorrected pixel-wise spectra (f).The fit of the adia SL spectrum was found to have a higher discrepancy between the fitted line and the measurement points compared to the other fitted spectra in this region.Substantial deviations between 0 and 1 ppm were found within the M T R asym measurements of adia SL in this ROI.Similarly, the M T R asym of the BPT revealed considerable deviations and a shift in the peak position.The boxplot highlights high asymmetric deviations around the median for adia SL and pronounced uncertainties for BPT saturation.display the Z-spectra for the two-pool simulation for a constant phase, while (e) and (f) display the spectra for an adjusted phase.The saturation of the BPT highly depends on the phase, for example, the z-magnetization at 1 ppm is 220 % higher than the simulation with the phase adjusted after each pause.These artifacts are periodically introduced, depending on the phase mismatch at the different off-resonances ∆ω making the spectrum non-smooth.Adjusting the phase after every pulse can avoid this, provided there are no field inhomogeneities.The OC pulse repeatedly reorients the magnetization vector in the z-direction after every block.This applies both to a constant phase and to a phase adjusted after each block, resulting in a nearly identical z-magnetization for the majority of off-resonant frequencies, which promotes a smooth spectrum (simulations of the magnetization vector over time for BPT and OC are displayed in the Supporting Information Figures S2,3).Analogously, in Figure 5 (g) and (h) after adding a B 0 shift of 0.1 ppm to the simulations with adjusted phase after every pause, artifacts are introduced into the spectrum of the BPT while the OC spectrum remains smooth under equal conditions.

Discussion
The aim of this work was to optimize the pulse shape of a CEST saturation pulse train to achieve high CEST contrast and stability in the presence of field inhomogeneities.The OC pulse train was compared against CW and Gaussian saturation in phantom measurements on a 7 'T preclinical scanner and against commonly used pulse types, such as Gaussian, BPT and advanced approaches like adiabatic SL pulses in phantom measurements on a 3 T clinical scanner.Additionally, Bloch-McConnell simulations were conducted to study the magnetization vector dynamics in the water pool during saturation with OC and BPT, to understand their behavior under the influence of B 0 inhomogeneities.The OC (Optimal Control) framework presented in this study successfully optimized the pulse shape.This enhancement not only improved contrast but also ensured stability amidst field inhomogeneities, outperforming other pulsed saturation strategies.
The different presented pulsed saturation regimes show substantial impact on the resulting spectra.
Since the presented pulse trains have the same B 1RM S , DC and saturation time, the only parameter that differentiates them is the time dependent B 1 (t) magnitude and phase.This work shows that with different pulse shapes several aspects of the CEST spectrum can be improved, suggesting the existence of an optimal pulse shape that enhances favorable features.Yoshimaru et al. were the first to successfully optimize pulse shapes, which were beyond the optimization of simple Gaussian or Sinc pulse shape parameters by applying a multi objective genetic algorithm [19].They compared the resulting pulse shapes against Gaussian and BPT / CW protocols in simulation and measurements.The primary focus of their research was on optimizing exchange weighting.Although they were successful in surpassing Gaussian, they were unable to approach saturation efficiency of CW/BPT saturation.
With the OC pulse designed in this work it was possible to outperform the exchange weighting of Gaussian pulse trains at 3 and 7 T in simulation and in creatine phantom measurements by up to 40-50 %.At 7T, the proposed OC pulse was able to reach the exchange weighting of CW saturation within a few percent and exceeded at 3 T BPT exchange weighting in all measurements.This performance difference also depends on the specific measurement parameters and physical environment at hand.In particular, the choice of the B 1RM S is of crucial importance as it influences the spillover and exchange weighting, which in principle depends on parameters such as the exchange rate, the relaxation times of water and CEST pool, B 0 and the resonance offset.In this study creatine was used for experimental validation.With an exchange rate of several hundred Hz and an off-resonance of roughly 2 ppm, it is a representative CEST agent in phantom studies and its physical properties lying within a typical range.The choice of B 1RM S value for our measurements was guided by the maximum M T R asym for CW saturation as stated by Zu et al. [18], a finding corroborated by our simulations (see Figure 1  We also conducted simulations to examine Fermi pulses, which, in terms of their shape, fall between Gaussian and block pulse trains (Supporting Information Figures S11 and S12).Compared to Gaussian saturation pulses, Fermi pulses demonstrated improved exchange weighting.However, their exchange weighting was notably less than that of CW and OC saturation.We observed that Fermi pulses introduced instabilities in the water peak within the -0.8 to 0.8 ppm range.These findings align with those reported in the study by Liu et al. [45].

Exchange weighting and stability in the presence of field inhomogeneities
The clinical utility of CEST effect measurements critically depends on their robustness to field inhomogeneities.Especially BPT saturation strategies can suffer from severe artifacts in the presence of B 0 inhomogeneities, which is especially obvious in the B 0 dependent signal behavior in Figure 3.
These artifacts could not be eliminated with B 0 correction and were only corrected to some extent by Lorentzian fitting.Adia SL saturation produced images with a notably higher contrast compared to Gaussian saturation, though it showed slightly less saturation than BPT.However, these images presented a subtle cloudy pattern, especially in the B 0 inhomogeneous region, a pattern that was further emphasized upon fitting.While adia SL saturation offers a valuable option for high-contrast and robust CEST imaging, OC saturation produced the most uniform images, coupled with the highest CEST contrast.
The robustness of the OC saturation to field inhomogeneities was evident in CEST creatine measurements at 3 T.In Figure 4 (c) the stability of the contrast was comparable for OC, adia SL and Gaussian saturation in the homogeneous B 0 region.The BPT saturation followed the B + 1 trend more closely indicating a stronger dependence to the B + 1 .However, as the line profile surpasses the center into a increasingly B 0 inhomogeneous region the robustness of the OC pulse became particular evident.While the contrast of the OC pulse was similar in the B 0 homogeneous and inhomogeneous region, the contrast generated by the other pulses showed deviations.The Gaussian profile was comparably robust, the adia SL saturation showed deviations with increasing B 0 shift and the BPT expressed sever artifacts starting at even slight inhomogeneities.This also became evident by investigating the statistical deviations in different ROIs (d).Over all ROIs the OC remained stable within all combinations of field variations.The susceptibility of the adia SL to strong B 0 inhomogeneities as described before by [21] became evident in ROI 2.Here the statistics as well as the spectra indicated that the adia SL saturation lead to an asymmetric spectrum around the water peak.This is especially challenging as artifacts in the M T R asym appear.Applying a Lorentzian fit to the water peak can yield inaccurate results, as a symmetric saturation profile is assumed.While the average spectra over the BPT were smooth, and the pixelwise fit was accurate, the deviation in the M T R asym and the statistics became apparent.The Gaussian saturation showed comparably high but in general stable deviations in all combinations of field inhomogeneities.
Remarkably, the OC pulses display inherent B 0 robustness, a feature not explicitly included in the cost functional but naturally emerging during the process, distinguishing our approach from previous studies like [24].Initially, the starting pulse, a BPT with constant phase, introduced severe artifacts/wiggles (Figure 5 b and c).Yet, the optimization found a solution that was effectively independent of phase changes, thus automatically compensating for unadjusted phase or B 0 offsets.
Incorporating an MT pool in the simulation similarly reduced the M T R asym for OC pulses, analogous to the effects observed with Gaussian and BPT saturation.

Insights from dynamic Bloch McConnell simulations
Dynamic Bloch McConnell simulations reveal that the origin of artifacts introduced by BPT saturation can be traced to a mismatch in phase between the pulses.After each block the magnetization vector ends with a non-zero magnetization in the x-and y direction.Only if the phase in the pause between the blocks is adjusted by ω n T pulse B 0 γ the resulting spectrum is smooth.However, B 0 inhomogeneities result in a phase mismatch after each block, leading to variable levels of saturation, depending on the off resonance.This causes artifacts in the spectrum (for more detailed explanation see Supporting Information Section 2).Analogously, similar artifacts can be observed in simulations by introducing a B 0 inhomogeneity of 0.1 ppm.These artifacts can be minimized by reducing the pause duration.However, the DC is restricted by SAR and hardware limits.Of particular interest is the ability of the OC pulse train to effectively flip the magnetization vector in the z direction, with M x and M y approaching zero after each pulse.Therefore, the subsequent pulse is independent of the previous pulse and independent of the pause duration.This makes the OC pulse inherently robust against B 0 inhomogeneities.
The reversal of the magnetization vector into the z direction after a saturation period is a concept specific to SL saturation techniques.Imaging with T 1ρ contrast instead of conventional CEST saturation has been shown to exhibit several advantages, such as high exchange weighting and selectivity [20,46,47].However, the off-resonant SL is susceptible to field inhomogeneities, underscoring the significance of Herz et al.'s recommendation to substitute tip-down and tip-up pulses with hyperbolic secant adiabatic pulses tailored to specific field strengths [21,22].Yet, the off-resonant SL saturation is elaborate to implement.At every off resonance ω n the phase of the tip up and tip down is different.Furthermore, the phase is also different if the frequency shift is either ω n > 0 ppm or ω n < 0 ppm.This is problematic if due to B 0 shifts a ω n > 0 ppm SL pulse is played at a ω n < 0 ppm point and vise versa.This makes the pulse vulnerable in the region around ω = 0 ppm, which leads to the asymmetric saturation as demonstrated in Figure 4.It appears that the other investigated pulses do not share this problem.Furthermore, to ensure robust tipping, the adiabatic pulses need to be long enough and they are generally SAR intensive.This limits the time and energy available for CEST saturation.
Interestingly, magnetization vector dynamics simulations suggest that OC pulses employ a saturation strategy similar to the off-resonant SL technique (this can be seen in the Supporting Information Figure S3).This process involves a tipping down phase, a saturation/locking phase, and a tipping back to M x , M y = 0 (for visualization see the videos of OC pulse saturation dynamic in the Supporting Information Video S1,2).Notably, OC pulses also achieve saturation during the tipping phase, thereby optimizing energy and time utilization.

Considerations for OC optimization
The simulations in Figure 1 indicate that the OC framework effectively generates pulse shapes for diverse parameter sets, producing results with exchange weighting and stability comparable to CW saturation.The pulse operates stable with changes in T 1 .The simulations indicate that the pulses can be used independently of the scanner field strength, which was demonstrated at a 3 T clinical system, and a 7 T preclinical system.This suggests that the optimized pulses are suitable for low to ultra high field strengths.
While the optimized pulses produce satisfying results for a broad range of parameters, they are not without limitations.Simulations indicate that high T 2 (> 250 ms) values and strong RF amplitude scaling (> 200 %) might introduce artifacts in the spectrum, although the M T R asym peak value remains stable.For applications with higher B 1 saturation, optimization of a pulse for a new target spectrum might be necessary.The wiggles at k = 2000s −1 in Figure 1 (c) can be removed by optimizing for a higher exchange rate (for details see Supporting Information Figure S4).
The pulses were optimized for a DC of 90 % with single pulse duration of 100 ms.As can be seen in simulation Figure 5 the pulse is robust against phase variations occurring during the pauses.This implies that the DC, and hence pauses, can be freely adjusted.However, a specific optimization is needed if a different pulse length t d is desired, as the single saturation duration is fixed here.
Chaining OC pulse trains can be used to increase total saturation time.Specifically, this allows for the total saturation time to be extended to integer multiples of the original T sat (for an example see Supporting Information Figure S10).
Optimization quality was improved by adopting real-valued pulses, which, unlike complex pulses, produce a symmetric spectrum, thereby reducing spectral artifacts.When given the freedom to optimize the phase, the optimizer found pulses that artificially enhanced the CEST effect.This means that we were able to measure a CEST peak in pure water even without a CEST agent.We are aware that this decision reduced the versatility of possible outcomes.Allowing phase optimization along the pulse train could potentially enhance desirable CEST saturation attributes, such as improved B + 1 robustness.Incorporating the phase into the optimization is a logical next step for further research.
The correct choice of the number of ∆ω in the z-spectrum was essential for the optimization.
Although the optimizer produced solutions fulfilling all requirements at the prescribed points, oscillations between the optimized ∆ω points could be observed with higher sampling rate.∆B 0 inhomogeneities displaced measurement points to areas with high oscillations, causing artifacts in the phantom spectra.However, these artifacts were minimized by employing a spectral resolution of at least 50 ∆ω points per ppm during optimization.

Conclusion
We have demonstrated that with the presented optimal control optimization RF pulse design method it is possible to design the pulse shape of an RF pulse for CEST imaging.On a 7 T preclinical scanner system, the optimized OC pulse train achieved saturation comparable to CW saturation and outperformed Gaussian saturation by up to 50 % in phantom measurements.The OC pulse train surpassed state-of-the-art CEST saturation strategies, producing the highest CEST

Figure 1
Figure 1 depicts the simulation results for the OC saturation for various parameter ranges and a comparison for M T R asym with Gaussian and CW saturation.Apart from the varied parameters the simulation parameters from the optimization were used.In (a) saturation was simulated for T 1 ∈ {0.2, 0.5, 1, 2, 4} s.The spectrum and the M T R asym peak (beneath) remained smooth and consistent for all T 1 times in the simulation.In (b) T 2 was varied for T 2 ∈ {0.01, 0.03, 0.08, 0.15, 0.25} s.

Figure 5 (
Figure 5 (a) depicts the BPT amplitude over time, while (b) shows the OC pulse train.(c) and (d) (e)).Recent in vivo measurements by Wang et al. 2024 for 3 T [44] further support this B 1RM S range.

1Figure 1 :
Figure 1: Simulation of OC Saturation for varying (a) T 1 times, (b) T 2 times, (c) exchange rates k, (d) field strengths B 0 , and (e) OC pulse amplitude scaling.The first row displays simulated spectra for OC saturation, while the second row shows the corresponding M T R asym .The third row presents the peak M T R asym values, comparing OC with CW and Gaussian simulations under the same physical conditions and pulse parameters.

Figure 2 :
Figure 2: Comparison of simulation and measurements in a creatine phantom at 7T. Simulated spectrum for CW, Gauss and OC saturation (a) with simulated M T R asym peak in (c).Uncorrected measurement spectrum for the same pulses (b) with corresponding M T R asym in (d).

Figure 3 :
Figure 3: For each pulse the (a) uncorrected M T R asym , (b) B 0 -corrected M T R asym , (c) TGVdenoised and B 0 -corrected M T R asym , (d) amplitude of the creatine CEST-peak in the two-pool-Lorentzian-fit and (e) TGV-denoised two-pool-Lorentzian-fit images is shown.An air bubble is visible at the top of the images introducing B 0 -inhomogeneities.(f) Vertical creatine CEST signal profile from top to bottom through the B 0 -inhomogeneities in (c).

Figure 4 :
Figure 4: OC, BPT, adia SL and Gauss creatine phantom CEST spectra measured at 3 T .(a) ∆ω shift map due to B 0 inhomogeneities with ROIs with different levels of field inhomogeneities.(b) Bloch Siegert B 1 -map divided by the nominal B 1 level.(c) Line profile in (b) through the center of the B 1 map and the most homogeneous B 0 region.(d) Boxplots of the the different ROIs in (a).(e) Mean spectrum in ROI 2. (f) Spectrum and fit of a pixel in ROI 2 (see arrow).Gauss) M T R asym of the spectra in (e).

Figure 5 :
Figure 5: Comparison of BPT (a) and the OC pulse train (b) simulations under different conditions: with constant phase (c,d), with an adjusted phase by ω n T pulse B 0 γ after every pause (e,f) and with an adjusted phase and an additional B 0 inhomogeneity of 0.1 ppm.
1,RM S = 1 µT.The B 1 level was selected based on simulations showing that both CW and Gaussian pulses approximately yield the maximum creatine CEST contrast for the given saturation time.The measurement spectra were compared to two-pool Bloch McConnell simulations.The parameters were the same as for the optimization, with the exception of the relaxation times, which were adjusted to the measured T 1 and T 2 times of

Table 1 :
Results from Figure 4 d, statistics for the median and IQR of the M T R asym in the ROIs 1-4 with different levels of B 0 inhomogeneities.