Parallel transmit pulse design for saturation homogeneity (PUSH) for magnetization transfer imaging at 7T

Purpose This work proposes a novel RF pulse design for parallel transmit (pTx) systems to obtain uniform saturation of semisolid magnetization for magnetization transfer (MT) contrast in the presence of transmit field B1+ inhomogeneities. The semisolid magnetization is usually modeled as being purely longitudinal, with the applied B1+ field saturating but not rotating its magnetization; thus, standard pTx pulse design methods do not apply. Theory and Methods Pulse design for saturation homogeneity (PUSH) optimizes pTx RF pulses by considering uniformity of root‐mean squared B1+, B1rms, which relates to the rate of semisolid saturation. Here we considered designs consisting of a small number of spatially non‐selective sub‐pulses optimized over either a single 2D plane or 3D. Simulations and in vivo experiments on a 7T Terra system with an 8‐TX Nova head coil in five subjects were carried out to study the homogenization of B1rms and of the MT contrast by acquiring MT ratio maps. Results Simulations and in vivo experiments showed up to six and two times more uniform B1rms compared to circular polarized (CP) mode for 2D and 3D optimizations, respectively. This translated into 4 and 1.25 times more uniform MT contrast, consistently for all subjects, where two sub‐pulses were enough for the implementation and coil used. Conclusion The proposed PUSH method obtains more uniform and higher MT contrast than CP mode within the same specific absorption rate (SAR) budget.


Introduction
Imaging at ultrahigh-field benefits from increased signal to noise ratio (SNR) 1 but is also hampered by larger transmit field ( 1 + ) inhomogeneity. This can lead to non-uniform excitation of the magnetization and undesired spatially varying contrast. Various hardware 2,3 and pulse design 4,5 solutions have been proposed and one of the most flexible is parallel transmit [6][7][8][9] (pTx) which uses multiple transmit channels to enable spatial and temporal manipulation of the 1 + field. Its most basic form, 'static' 1 + shimming [10][11][12] , attempts to create a more homogeneous 1 + field by applying amplitude and phase weightings to the individual transmit channels without changing the radiofrequency (RF) pulse waveforms themselves. Much greater control can be achieved if we instead consider the magnetization rotation using both RF and gradients; in this case a desired 'flip angle' distribution is usually designed by tailoring the RF pulse waveforms on each channel. This is often achieved using the small tip angle approximation in which case the design problem can be considered using the excitation k-space concept 13 . For larger rotations this concept breaks down but various methods still exist to design RF pulses by considering the rotation of magnetization directly 14- 16 .
Although a great diversity of RF pulse design methods exist, they typically have in common the use of the Bloch equation 17 to model the effect of the applied RF and/or gradient fields. Assuming short duration pulses, the effect of applying these fields is to rotate the magnetization. The Bloch equation can successfully model the magnetization dynamics of free water, but in biological tissues there is usually also a significant pool of semisolid magnetization 18,19 . This latter pool is affected differently by RF fields and can exchange magnetization with the free water, a phenomenon usually referred to as Magnetization Transfer 20,21 (MT). Two common assumptions of the semisolid pool are: (1) it can be modeled as having no transverse magnetization due to its very short 2 ≈ 10 ; (2) its longitudinal magnetization directly saturates at a rate proportional to the applied RF power 22 (i.e. | 1 + | 2 ). These properties and the coupling with free water magnetization can be modeled using the so-called binaryspin-bath (BSB) model 20 .
Since the saturation of semisolid magnetization depends on | 1 + | 2 and not simply 1 + , the effect of 1 + inhomogeneities is more severe. Existing RF pulse design methods that might be used to correct for non-uniform flip angles in free water magnetization will fail if used for designing semisolid saturation pulses because: i) the semisolid has no transverse component that can be rotated by any applied gradients that are often used to improve excitation properties 23,24 and ii) the saturation rate of its longitudinal magnetization depends on | 1 + | 2 and not 1 +20, 22 .
An important distinction here must be drawn between the saturation of semisolids, which is the subject of this work, and general 'saturation' pulses that are used (often with spoiler gradients) to suppress magnetization from free water and/or solutes. For the latter type there are examples using standard pTx pulse design methods 25,26 , as there is transverse magnetization amenable to rotation from the RF pulses and gradients. For the remainder of this article the term 'saturation' is used to refer to semisolid saturation, unless otherwise specified.
In this work we propose a novel RF pulse design framework for semisolid saturation, called PUlse design for Saturation Homogeneity (PUSH). We first explore a general case of RF pulse design in the presence of semisolids, and then propose a simple exemplar method using trains of short sub-pulses with pTx and demonstrate the efficacy of this approach for MT-weighted imaging at 7T.
On the other hand, systems with MT can be described by the BSB model 22 which contains two pools corresponding to free water ( ) and semisolid ( ) magnetization: The operators ̃, ̃ and ̃ include the effects of RF, gradients, relaxation, and exchange. In these expressions 1 (=1/ 1 ) is the semisolid longitudinal relaxation rate, M 0 is the semisolid equilibrium magnetization, is the exchange rate from M to M (vice-versa for ), and 〈W〉 is the average saturation rate 22 that models the semisolid response to RF. In case of RF irradiation at a single offresonance frequency then 〈W〉 is given by 22 where 〈| 1 + | 2 〉 is the mean squared 1 + over pulse duration , and is the semisolid absorption lineshape that depends on its transverse relaxation time 2 and on the frequency shift − γΔB z . Typically the absorption lineshape has much broader bandwidth 21 than the RF or 0 field variations in the absence of gradients, such that ( − γΔB z , 2 ) ≈ ( , 2 ). It has also been observed that may have a chemical shift away from water (e.g., Jiang et al 27 observed ≈ −2.6ppm in white matter)this shift should be considered part of the definition of . Although Eq. [3] defines 〈W〉 for single frequency irradiation, it can also be calculated in some cases for RF pulses with multiple frequencies (e.g., multiband pulses 28 ).

RF pulse design
In solving the Bloch equation (Eq. [1]) for a short RF pulse, matrix and vector can be neglected as relaxation typically occurs over a longer timescale, hence the magnetization dynamics comprises of rotations determined by . This can be solved by discretizing the sequence parameters in constant piecewise timesteps of duration Δ , each producing a rotation ( ): ( + Δ ) = exp( ( )Δ ) ( ) = ( ) ( ) [4] whereas the full rotation full of the magnetization can be calculated by taking left-wise multiplication over all ( ).
Similarly, in the BSB model (Eq. [2]) relaxation and exchange can be assumed to occur over a longer timescale than the typical RF pulse allowing them to be neglected for RF pulse design. The magnetization response is thus determined by ̃ and can also be solved by discretizing time: In the absence of exchange the response of the free water and semisolid pools is decoupled (̃ is block diagonal), thus the matrix exponential is the exponential of its diagonal terms. The full magnetization response to an RF pulse can then be calculated by taking the product of the matrix exponentials from all timesteps over the duration of the RF: where the free water and semisolid pools responses are independent from each other. Design of RF pulses for free water magnetization usually target a desired flip angle from rotation matrix full . On the other hand, to control saturation of semisolid magnetization we can target a desired average saturation rate 〈W〉 . Thus, the RF pulse design for both pools can be cast as a joint optimization: where ∈ [0,1] balances the error between the two terms and and are the RF and gradient waveforms, respectively.

Pulse design for Saturation Homogeneity
In this work we explored design of RF pulses to achieve uniform semisolid saturation, using pTx systems. To this end we have considered single frequency high power saturation pulses 29,30 ( Figure 1A) applied at large offset frequency . We assume that is much larger than the saturation pulse bandwidth such that we can neglect the effect on the free water magnetization, i.e., it has null flip angle ( ≈ 0)effectively performing the pulse design in Eq. [7] with = 1. Furthermore, according to Eq. [3] for single frequency irradiation we can control semisolid saturation using 〈| 1 + | 2 〉 instead of 〈W〉 in the pulse design. For a pTx system the applied 1 + field is the linear superposition of the fields from its ℎ channels: [8] where ( ) are the transmit sensitivity maps (units of /V) and ( ) are the RF waveforms for each channel (units of V). In this implementation we employed short repetition time (TR) sequences in which saturation depends on the cumulative effect over many TR periods, scaling with the mean squared 1 + (Eq. [3]) averaged over the TR 31,32 instead of the pulse duration : Here the contribution from any other pulses during the same TR period (e.g., excitation pulse as described in Methods) is neglected as they typically have much less power than the designed saturation pulse. The advantage of designing for 1 (= √〈| 1 + | 2 〉) at the sequence rather than pulse level is that the former is typically the limiting factor for an MT-weighted sequence since it scales with the SAR; exposing this limit allows more flexibility to optimize the sequence within this constraint.
In this work we applied the same normalized waveform ( ) (arbitrary units) in each channel scaled by a complex weight w (units of V), moving the sum outside the integral: where 〈 2 〉 is the mean squared 1 + of the normalized waveform. Similarly to spokes 33,34 /kT-points 35 , the pulse can be extended by concatenating sub-pulses, designated as PUSH-, with each subpulse weighted differently: where is the sub-pulse index. Each sub-pulse produces its own spatial mean squared 1 + , such that the total 〈| 1 + | 2 〉 is the sum of the contributions from all sub-pulses (example in Supporting Information Figure S1). Finally, the RF complex weights w can be designed to achieve a desired saturation by solving the optimization: where β( ) specifies the desired spatial √〈| 1 + | 2 〉 = √〈W〉/ 2 ( , 2 ). Note that the optimization has been rewritten in terms of the square root of 〈| 1 + | 2 〉doing so does not change the global optima.
This was done so that the special case of a saturation pulse consisting of a single sub-pulse would reduce to Magnitude Least-Squares (MLS) 1 + shimming 11,12,36 of the saturation pulse. The terms 'PUSH-1 (1 sub-pulse)' and 'static shimming' will be used interchangeably from here onwards. Optimization was constrained to be within local SAR limits for a total of N VOP VOPs 37 , as well as average power per channel and maximum voltage V max per sub-pulse and per channel 38 . In the case of using the circular polarized (CP) mode, the complex weights are defined as w = w CP exp(− 2 ( − 1) ℎ ⁄ ), where w CP was determined by minimizing the cost function (Eq. [12]) and 2 = −1 denotes the imaginary unit.

Methods
All experiments were performed using a 7T scanner (MAGNETOM Terra, Siemens Healthcare, Erlangen, Germany) in prototype research configuration, with an 8Tx/32Rx head coil (Nova Medical, Wilmington MA, USA).

Pulse sequence setup
To illustrate the PUSH concept, we used a simple MT-weighted spoiled gradient-recalled echo (SPGR) sequence containing one saturation and one excitation pulse per TR ( Figure 1A). The saturation subpulses ( Figure 1B) applied at offset frequency = 2kHz had a Gaussian waveform (Time Bandwidth Product= 2.27, = 4 ) and its complex weights were determined using Eq. [12]. Pulse optimization was solved in Matlab (Mathworks Inc., Natick, MA) using the interior-point algorithm from the fmincon routine, providing first and second order derivatives; constraints (SAR 10g,max = 20W/kg over 8 VOPs provided by the vendor and in first level SAR mode 39 , V max = 207V, max = 24W) were evaluated within the vendor pulse design framework included in the scanner console software (release Syngo.MR VE12U). A multi-start strategy with 10 random seeds proved to obtain consistent solutions. As we focused on the saturation pulse design, the excitation was always in CP mode. To minimize the impact of the excitation pulse on the MT contrast, the flip angle was minimized balancing SNR and its inhomogeneity profile (Supporting Information Figure S2). This way the excitation pulse also had negligible power compared to the saturation pulse such that √〈| 1 + | 2 〉 in Eq. [12] can be assumed equivalent to the sequence 1 .

Simulations
To explore the pulse design performance, saturation pulses with different number of sub-pulses (1, 2 and 3) were designed offline for a spatially invariant β ranging from 0.1 to 2 in steps of 0.1 . The optimizations were performed for both 2D axial slices and 3D volume of brain transmit maps from an 8-channel pTx system (details below), with each 2D single slice/3D optimization taking ≈ 7/22 seconds in Matlab R2020b (Mathworks Inc., Natick, MA) on a desktop computer (Intel i9-10900X @ 3.70GHz, 64GB of RAM, not parallelized). The solutions were analyzed in terms of their 1 maps and normalized root mean square error (NRMSE).
To predict the impact of spatial variation in mean squared 1 + on the MT contrast, Magnetization Transfer Ratio (MTR) maps were simulated using the definition: where and are the steady-state signals acquired with and without the saturation pulse, respectively. For the MTR simulations the steady state of an SPGR sequence was calculated by solving Eq. [2] assuming the whole brain to have uniform tissue parameters similar to white matter 40 1.85s −1 , 2 = 9.6μs and a Super-Lorentzian absorption lineshape (centered at -773Hz) 27 . Different saturation pulses optimized offline for β = 1 were applied combined with a small excitation flip angle of = 5 ∘ .

Experiments
In vivo scanning of five healthy volunteers was performed in accordance with local ethical approval. MTR maps were acquired for 2D and 3D imaging, as described in the subsections below. The saturation pulse was designed online as described in subsection 2.3 with calculation fully scanner-integrated within a Matlab R2012b (Mathworks Inc., Natick, MA) framework from the vendor, taking ≈ 30 seconds. Prior to pulse design off-resonance Δ 0 mapping was performed using a dual-echo SPGR sequence and 1 + mapping was performed using a pre-saturation turbo-FLASH sequence 42 (both with resolution 4 × 4 × 6 mm 3 ); 1 + was corrected for bias using an empirically determined correction factor 43 . The reference voltage (V ref ) determined by the scanner's built-in adjustment steps was also recorded: a higher reference V ref indicates lower efficiency in generating 1 + , and hence lower achieved 1 + for a given SAR level. A signal intensity-based mask from the vendor's framework was used to not impair the workflow but was pre-processed to remove non-brain tissue voxels by eroding each axial slice with a 3-pixel (12mm) radius disk and cropping axial slices that included voxels from the mouth and jaw. Prior to MTR maps calculation (Eq. [13]), images were registered together using FSL BET 44 and FSL FLIRT 45 , and white matter (WM) segmentation was performed using FSL FAST 46 , further eroded with a 1-pixel radius disk to reduce partial volume effects. 1 maps were simulated retrospectively using the 1 + maps and the pulses optimized online (during the scan).

2D imaging
2D MTR maps were acquired in all subjects for a single axial slice in the middle of the brain (resolution 1 × 1 × 5 mm 3 , matrix size 220 × 220, TR=22ms, TE=4ms, BW = 220Hz/Px, 4 averages). Data were acquired using three different saturation pulses (CP mode, PUSH-1, PUSH-2) and for four β (0.7 , 1.0 , 1.3 , 1.6 ). For each MTR map a set of and images were acquired, firstly with 30 seconds of dummy pulses (to stabilize the RF output) followed immediately by with 10 seconds of dummy pulses, resulting in = 1: 19s per MTR map. Fewer dummy cycles were required for since it was acquired immediately after . White matter segmentation for further analysis was performed using the MTR map obtained with PUSH-2 at β = 1.3 due to its uniform contrast (as shown later).
An additional MP2RAGE 48 acquired at the same resolution was used for segmentation.

Gradient blip experiment
To explore the impact of applying gradients between RF sub-pulses on the semisolid saturation, we carried an experiment adding gradient blips to the PUSH saturation pulse trains. These gradients would affect the 'flip angle' if this pulse were applied to free water magnetization, but not the semisolid saturation if it behaves as modeled with longitudinal magnetization only. 2D MTR maps (same protocol as in 3.3.1 with TR=27ms) were acquired in one subject with saturation pulses designed using PUSH-3 for β = 1 , both excluding and including gradient blips between the sub-pulses. For the latter the gradient blips were 100 s in duration in the x and y-directions, each producing a 4 phase roll across the FOV. The 1 and flip angle maps of each pulse were computed (flip angle calculated at the RF offset frequency with Eq. [6]) and compared to the measured MTR maps. Figure 2 examines the error in 1 as a function of the target β for 2D optimization on a middle axial slice ( Figure 2A) and 3D whole-brain optimization ( Figure 2B). Supporting Information Figure S3 further expands on Figure 2A for other axial slices. In all cases CP mode has a constant NRMSE until reaching the local SAR limits when the error increases as the voltage is capped. On the other hand, PUSH-1, -2 and -3 perform better than CP mode across all β for both 2D and 3D imaging. For the 2D case PUSH-1 (i.e. static shimming) gives an NRMSE 2 times smaller than CP mode for β ≤ 0.4μT in the middle slice ( Figure 2A) but its performance worsens as β increases. PUSH-2 and -3 perform equally well, with a constant NRMSE 6 times smaller than CP mode (for β ≤ 1 ) in the middle slice.

Simulations
Remarkably, PUSH-2 and -3 still perform well for β ≥ 1μT, beyond where CP mode reached the local SAR limits. Generally, in 2D the performance gain offered by PUSH was larger for middle and inferior slices while superior slices had less inhomogeneity to begin with. For 3D imaging ( Figure 2B) similar trends were observed, but with overall larger NRMSE. For β ≤ 1 PUSH-1 obtains an NRMSE ≈ 17% smaller than CP mode, whereas PUSH-2 and -3 perform again equally well and still achieve an NRMSE 2 times smaller than CP mode. Figures 3 and 4 show computed 1 maps for 2D and 3D imaging, respectively, for some of the solutions in Figure 2. Note that in Figure 4 some stripes are visible in the sagittal and coronal planes; these were found to be caused by artefacts in the acquired 1 + maps. In Figure 3 the 1 for CP mode scales up with β and caps after reaching the local SAR limits. PUSH-1 achieves more uniform 1 for β ≤ 0.7μT but gets progressively less homogeneous with increasing β, producing solutions with "holes". PUSH-2 and -3 achieve 1 maps similar to one another that are more uniform up to larger β. For 3D imaging, Figure 4 shows that CP mode produces a 1 pattern with center brightening. The 1 produced by PUSH-1 is very similar to CP mode with slight improvements for β ≤ 0.7μT. With PUSH-2 and -3 the 1 is more uniform in the middle slices and up to larger β, however both solutions underdeliver in the superior and inferior slices of the brain.  1 is shown below the respective maps.  Figure 5 shows simulated 2D and 3D MTR maps for β=1 solutions from Figure 2. Note that all simulated MTR maps have been calculated using white matter properties over the whole brain: they are intended to visualize the spatial variations of saturation rather than the actual MTR contrast that will be seen in a scan. Both CP mode and PUSH-1 solutions yield similar MTR maps with center brightening for both 2D and 3D. However, for 2D imaging PUSH-1 can also yield solutions with contrast "holes" for some slices. On the other hand, PUSH-2 and -3 achieve similar strong improvement in 2D and 3D, with a standard deviation ≈ 2 times smaller than CP mode. A drop in achieved MTR towards the superior and inferior slices is seen in 3D. For all cases the MTR maps have a good correlation with the respective 1 in Figures 3 and 4. Figure  2 for =1 as the saturation pulses. For 2D slices 6, 10, 12 (slice in Figure 3), 14 and 18 (from left to right) that were individually optimized are shown. The green arrow in (A) points to a contrast "hole" seen in some slices optimized with PUSH-1. The mean ± standard deviation of MTR over the whole volume is shown below the respective maps.

2D imaging
2D MTR maps from one subject are shown in Figure 6 for different pulses and β. CP mode shows a constant contrast pattern brighter in the center, that scales up with the β until it reaches the local SAR limits, after which the voltage is capped and no more RF power is delivered with increasing β. PUSH-1 achieves more uniform contrast for the smallest β but then resembles CP mode up to the largest β where it has a "hole" in the contrast whilst increasing the contrast everywhere else. PUSH-2 yields uniform contrast for all β with up to 4 times smaller dispersion except for β = 1.6 , where it is slightly less bright in the center. The MTR maps correlate very well with the corresponding 1 maps in Supporting Information Figure S4.   Figure 7 illustrates the MTR distribution in white matter for 5 subjects as a function of the saturation pulse and β; the subjects are arranged in order of increasing V ref . The MTR distributions are consistent across all subjects, with PUSH-2 yielding narrower distributions for all β. PUSH-1 has narrower distributions for the smallest β, but at 1.6 its distribution exhibits a heavy tail towards low MTR values as its maps have "holes" (Figure 6). Nevertheless, both PUSH-1 and -2 can achieve higher mean MTR than CP mode for the largest β, as in all subjects CP mode reached the local SAR limits below 1.3 . The MTR in Figure 7 correlates well with the corresponding 1 distributions in Supporting Information Figure S5. The subjects are ordered by increasing V ref , showing an expected negative correlation between MTR and V ref for the largest β (also illustrated by Supporting Information Figure  S6), as higher V ref means lower maximum 1 + peak.    Figure 9 shows results from PUSH-3 pulses without ( Figure 9A) and with ( Figure 9E) gradient blips between sub-pulses; MTR maps were acquired while 1 and flip angle maps were simulated from acquired 1 + and Δ 0 maps. The MTR maps ( Figure 9D and 9H) are virtually identical and very uniform, in agreement with their respective 1 maps (Figures 9B and 9F). The flip angle maps ( Figures 9C and 9G) computed by also considering rotation induced by the gradients have a very different appearance and are both non-uniform.

Discussion
This work presents a novel pTx pulse design to overcome 1 + inhomogeneity in MT imaging at ultrahigh field by controlling semisolid saturation through the mean squared 1 + . This can be performed either instead of or in addition to controlling the excitation properties of water for a given RF pulse; in this work we focused on testing the former case which is relevant to application of off-resonance saturation pulses. PUSH was tested in simulations and in vivo, yielding more uniform MT contrast. Current pTx pulse design methods 24,33-35 usually employ a combination of RF and 0 gradients to optimize the flip angle of the resulting excitation. These methods are unsuitable for designing semisolid saturation pulses for MT imaging since the semisolid pool has no transverse magnetization and instead saturates directly with | 1 + | 2 . This was experimentally confirmed (Figure 9) by applying the same RF pulse twice, with and without gradient blips in-between its sub-pulses. The gradients blips do not change | 1 + | 2 but drastically alter the flip angle if applied to free water magnetization. The measured MTR maps show no difference in the MT contrast, supporting the fact that 'flip angle' is not a useful metric to use when designing or describing semisolid saturation pulses. A simpler alternative to 'dynamic pTx' pulse design is 1 + shimming 11,12,36 which aims to create a spatially uniform 1 + distribution. An optimal 1 + shimming solution would also achieve a uniform | 1 + | 2 meaning that in principle 1 + shimming is a special case of PUSH where the RF pulse is 'static' (i.e. pTx degrees of freedom are not modulated through the pulse). To connect PUSH with 1 + shimming, the optimization (Eq. [12]) is formulated in terms of the square-root of 〈| 1 + | 2 〉, such that with 1 sub-pulse it simplifies to a magnitude least-squares 11,12,36 1 + shimming design.
A related alternative approach is the MIMOSA 49 pTx design proposed to homogenize saturation in pulsed CEST. In that case a train of saturation pulses interleaved with spoiling gradients can be approximated by an equivalent continuous wave saturation whose effective 1 + is better described by the 1 50 over the pulse train. Hence the MIMOSA design applies two complementary modes 51,52 with the objective of homogenizing 1 for CEST saturation. This effectively results in a solution equivalent to PUSH-2, but for the special case of using two pre-defined modes of the pTx coil.

PUSH performance for MTR imaging at 7T
Simulations show that in both 2D and 3D imaging (Figure 2) PUSH-2 and -3 yield a strong improvement in the homogeneity of 1 for all β achievable (within SAR limits) with CP mode, whereas with PUSH-1 (static shimming) the improvements reduce as β increases. Remarkably, PUSH-2 and -3 sustain these improvements beyond β achievable with CP mode, delivering higher and more uniform 1 . For 2D imaging (Figure 3 and Supporting Information Figure S3) these improvements are more substantial in the inferior and middle axial slices of the brain, however the maximum β achievable with CP mode is considerably smaller for the inferior slices. For 3D imaging ( Figure 4) these improvements are smaller, with 1 decreasing from the middle towards the inferior and superior slices. In all simulations 2 or 3 sub-pulses yield very similar results, suggesting that 2 sub-pulses are enough to explore the variability in the transmit sensitivity maps used, hence in most in vivo experiments a maximum of 2 sub-pulses was used.
The 2D in vivo experiments shows up to 4 times more uniform MTR maps with PUSH-2 ( Figure  6), corroborated by the narrow distributions of MTR in WM (Figure 7). Higher maximum MTR is obtained with PUSH-2, as expected from the simulations. Moreover, PUSH-1 (static shimming) solutions for β = 1.6 have pathological contrast with "holes", which is known to affect shimming solutions 53 . The 2D results are consistent across all subjects, with the MTR maps correlating very well with the respective 1 maps. Some inter-subject variability is observed for the maximum MTR achieved ( Figure 7D), relating to the reference voltage of each subject (Supporting Information Figure  S6). This is understandable since a higher reference voltage indicates the subject experiences a higher SAR per unit of achieved 1 + leaving less room for optimization.
The 3D in vivo experiments also show more homogeneous MTR maps with PUSH-2 ( Figure 8) but with a more modest 25% improvement in homogeneity. The distribution in WM shows smaller MTR values in the top and bottom slabs, agreeing with the respective 1 distribution (Supporting Information Figure S7). This effect is also observed in the 3D simulations ( Figure 4). A potential halfway point between the 2D and 3D results would be to use a multi-slab approach where saturation pulses are designed separately for each slab (though they are spatially non-selective due to the broad semisolid lineshape) and are paired with slab selective excitation pulses.

Impact of RF coil design
While current pTx methods can employ gradients to enhance spatial encoding of RF pulses designed to achieve rotations of magnetization, PUSH relies solely on the transmit sensitivity maps to homogenize the MT contrast. This is seen particularly in the performance of PUSH in 3D where there is a persistent decrease in the achieved MTR in the superior and inferior regions. This is consistent with the limited coverage and lack of pTx control over 1 + variation in the z-axis (head-foot) from the circumferentially arranged transmit elements in the coil used. It is likely that the proposed method would benefit from alternative coil geometries 10 , e.g. more channels and/or different distribution, to achieve a greater control of the mean squared 1 + spatial distribution.
For the coil used in this work we found that more than 2 sub-pulses does not improve the MT contrast homogeneity, but this may also prove not to be the case for alternative coil designs. More subpulses might also be beneficial in the case where peak voltage is the active constraint, whereas in the current implementation with the sequence and hardware used, local SAR was always the most limiting.

Assumptions and future extensions
Although exchange (Eq. [2]) is neglected over the RF duration, its cumulative effect over the whole pulse sequence makes MTR sensitive to both saturation of the semisolid and rotation of the free water 54,55 . Thus the excitation pulse can also affect MT contrast because i) it applies some power and ii) it rotates the free water magnetization that exchanges with the semisolid pool. In order to focus only on the saturation pulses our experiments employed a low excitation flip angle to minimize MTR sensitivity to excitation inhomogeneities, as excitation pulses used CP mode. As a result, the observed MTR is highly correlated with the 1 over the entire sequence; the contribution of excitation pulses to the 1 is negligible. This simple embodiment is used as a means to illustrate the key concept; however a future implementation might also consider designing uniform excitation pulses using methods such as kT-points 35 or SPINS 24 potentially as part of a joint optimization problem (Eq. [7]). Likewise, it is not necessary to compute pulses in terms of the sequence 1 as done here. Use of the root-mean squared instead mean squared 1 + in Eq. [12] has the advantage that in the case of 1 subpulse it simplifies to a static shimming problem, making it a special case and allowing for a direct comparison. Likewise, averaging over the sequence TR rather than the pulse duration connects more closely to the expected SAR limits 31 , and MT contrast for sequences with short TR where the continuous wave approximation is still valid 31,32 . However, in sequences with long TR this approximation breaks down and different exchange times affect MT contrast, so it is more appropriate to consider 〈| 1 + | 2 〉 over the pulse duration (as in Eq. [3]) instead of over the TR (Eq. [9]). Although gradient blips are observed not to affect the MT contrast, gradients applied during (as opposed to in between) RF pulses are expected to affect the semisolid saturation. According to Eq. [3] the saturation depends on 〈| 1 + | 2 〉 but also on the absorption lineshape ( − γΔB z , 2 ). Thus, theoretically it is also possible to control the saturation using applied gradients, which could be an avenue to explore, though this would require prior knowledge of the absorption lineshape 21,56 .

Conclusion
This work proposed a novel RF pulse framework called PUlse design for Saturation Homogeneity (PUSH) for design of RF pulses considering their saturation effect on semisolid magnetization relevant to magnetization transfer imaging. It was also demonstrated that adding gradient blips between RF sub pulses as commonly used by standard pTx methods does not affect the MT contrast; the 'flip angle' of a saturation pulse is not a meaningful way of describing its operation.
The specific case demonstrated in this work was the design of off-resonance saturation pulses where on-resonance effects can be neglected. Simulations and in vivo experiments showed that for the 8 channel RF coil used in this work PUSH can obtain up to 4 and 1.25 times more uniform MT contrast in 2D and 3D imaging, respectively, achieving monomodal distributions of MTR that correlate very well with the corresponding applied 1 . Moreover, PUSH delivered higher 1 than CP mode under the same SAR budget, thus also obtaining stronger contrast. Figure S3: (A) NRMSE of B rms 1 for the axial slice positioned as indicated by the red line in the sagittal plane in (B), comparing CP mode with the optimized PUSH solutions using 1, 2 and 3 sub-pulses (curves for 2 and 3 sub-pulses are superimposed due to nearly identical performance). The gray area represents β where CP mode reached the local SAR limits and its voltage is capped. Slice 12 corresponds to the solution in Figure 2A. To navigate through different slices this document needs to be open on a JavaScript-supporting PDF viewer, such as Adobe Acrobat Reader.

Supporting Information
Supporting Information Figure S4: Corresponding 2D B rms 1 maps for the MTR maps in Figure 6. Different rows correspond to different pulses (top row: CP mode; middle row: PUSH-1; bottom row: PUSH-2) and columns correspond to different β (increasing from left to right). Below each is the mean ± standard deviation of B rms