Steady‐state imaging with inhomogeneous magnetization transfer contrast using multiband radiofrequency pulses

Purpose Inhomogeneous magnetization transfer (ihMT) is an emerging form of MRI contrast that may offer high specificity for myelinated tissue. Existing ihMT and pulsed MT sequences often use separate radiofrequency pulses for saturation and signal excitation. This study investigates the use of nonselective multiband radiofrequency pulses for simultaneous off‐resonance saturation and on‐resonance excitation specifically for generation of ihMT contrast within rapid steady‐state pulse sequences. Theory and Methods A matrix‐based signal modeling approach was developed and applied for both balanced steady state free precession and spoiled gradient echo sequences, accounting specifically for multiband pulses. Phantom experiments were performed using a combination of balanced steady state free precession and spoiled gradient echo sequences, and compared with model fits. A human brain imaging exam was performed using balanced steady state free precession sequences to demonstrate the achieved contrast. Results A simple signal model derived assuming instantaneous radiofrequency pulses was shown to agree well with full integration of the governing equations and provided fits to phantom data for materials with strong ihMT contrast (PL161 root mean square error = 0.9%, and hair conditioner root mean square error = 2.4%). In vivo ihMT ratio images showed the expected white matter contrast that has been seen by other ihMT investigations, and the observed ihMT ratios corresponded well with predictions. Conclusions ihMT contrast can be generated by integrating multiband radiofrequency pulses directly into both spoiled gradient echo and balanced steady state free precession sequences, and the presented signal modeling approach can be used to understand the acquired signals.


| INTRODUCTION
Inhomogeneous magnetization transfer (ihMT) is a method for generating contrast specific to substances with lamellar or otherwise ordered microstructure that can support local nonzero dipolar magnetization order. [1][2][3] It has been demonstrated that this effect could have high specificity for myelinated tissue, and ihMT methods are now being established as structural markers for myelin 4 with potential use for studying demyelinating conditions, such as multiple sclerosis. 5 The ihMT effect results in different observable MT characteristics when single or dual-frequency off-resonance irradiation are used during imaging. Novel pulse sequences incorporating either multiband or rapidly alternating singlefrequency saturation pulses have been proposed to generate ihMT weighted images 6,7 demonstrating strong white-matter specificity. More recently, quantitative approaches to measure the dipolar relaxation time 8 and influence of ihMT on free water T 1 have also emerged. Existing work has focused on a preparation-based paradigm in which saturation pulses are interleaved with readout periods for measurement, as illustrated by Figure 1. Classic pulsed MT sequences 9 often consist of separate saturation and excitation pulses preceding a gradient echo readout ( Figure 1A), while a recently proposed ihMT sequence 7 uses multiple saturation pulses with fixed or alternating frequency offsets followed by multiple gradient echoes ( Figure 1B). In our own previous work, 10 we have developed nonselective multiband excitation pulses to control for MT effects in variable flip angle relaxometry, demonstrating that use of constant radiofrequency (RF) power over all flip angles leads to more stable relaxation time measurements. This sequence, illustrated in Figure 1C, is essentially a modified pulsed MT sequence with nonselective multiband pulses used to simultaneously perform saturation and excitation, and may be combined with a spoiled or balanced readout. Multiband pulses potentially lead to shorter repetition times (TRs) and include the flexibility to add either single or dual frequency saturation bands for generation of ihMT contrast.
In this work, we explore the use of the proposed multiband sequence for generation of ihMT contrast and formulate a matrix-based signal model for ihMT that may be used to directly calculate steady-state signals for spoiled gradient echo (SPGR) and balanced steady state free precession (bSSFP) sequences when dipolar terms are present and multiband F I G U R E 1 A, Classic pulsed MT weighted sequence containing off-resonant saturation pulse followed by gradient echo readout. Here shown as a 3D sequence, although this can also be made slice selective. B, Pulsed ihMT sequence from Mchinda et al 7 consisting of multiple saturation pulses using either fixed or alternating frequency offsets for generation of ihMT contrast, followed by multiple gradient echo readouts. C, Proposed sequence using nonselective multiband pulses that simultaneously saturate semisolid magnetization components while exciting water on-resonance. This may be combined with either balanced (illustrated) or spoiled readouts. Each of these diagrams shows 1 TR period of a steady-state sequence pulses are used. We investigate the validity of this matrixbased model, compare predictions with phantom experiments, and present in vivo ihMT ratio images formed using the proposed sequences.

| Bloch-McConnell-Provotorov equations
The commonly used "binary spin bath" model for magnetization transfer effects divides tissue magnetization between free water and semisolids, labeled with superscripts f and s, respectively. In this model, the ( whose temporal evolution is governed by: Matrices and C describe evolution in the absence of RF: where R f 1,2 are the free water relaxation rates (inverse of relaxation times); R s 1Z and R s 1D are semisolid Zeeman and dipolar relaxation rates, respectively; Δ is the off-resonance frequency; k is the exchange rate constant between free water and semisolid pools; and M f 0 and M s 0 are thermal equilibrium magnetizations for the free and semisolid pools, respectively. Note that, in this article, k represents the intrinsic exchange rate as opposed to a directional rate (i.e., forward or reverse) as sometimes used in the literature (for example, see Yarnykh). 13 For simplicity, we assume that k and R s 1Z are the same for pools s1 and s2, and use this formulation in this article. A more general version is given in the Appendix. Pools s1 and s2 occupy fractions 1 − and , respectively, of the total semisolid compartment M s 0 and we adopt the convention M f 0 + M s 0 = 1. Figure 2 illustrates the relationships between compartments in this model.
consists of free water in contact with a semisolid compartment that has two subparts s1 and s2. The former is a "classic" MT model with only Zeeman ordered magnetization; the latter contains a dipolar ordered component, which only exchanges with the Zeeman order under the action of off-resonant RF saturation To handle multiband pulses, we explicitly decompose them into contributions from separate frequency bands: where Δ j is the (angular) frequency offset of the j th band and b{Δ j ; t} is the time domain pulse waveform associated with component frequency Δ j . It then follows that matrix describing interactions due to RF pulses may be written as: free is simply the Bloch equation for the free pool. For the multiband pulses used in this work, only the onresonance component b{0; t} need be considered. semi describes absorption by the semisolid pools and is explicitly written as a sum over frequency bands as suggested by Goldman 11 for the case of "fast modulation". The RF saturation rate is defined as where g Δ, T s 2 is the absorption lineshape of the semisolid and T s 2 is its transverse relaxation time. The term Δ � j ≡ Δ j −Δ accounts for off-resonance effects on the local absorption response; for simplicity, these have been neglected in this work and we have instead assumed that Δ � j ≈ Δ j . The term loc characterizes the strength of local field fluctuations, which relate to dipolar interactions; this may be derived from the second moment of the lineshape. 1,14 As with exchange and relaxation rates, we assume both semisolid compartments have the same absorption lineshape.
There are a few properties to note for sequence modeling: (i) the first element of semi is exactly as commonly used to describe "classic" MT effects; (ii) exchange between M s2 Z and M s2 D is purely RF driven, while exchange between M f z and M s1,2 Z occurs continuously; (iii) coupling terms in semi are linear in Δ j so equal off-resonance RF power applied simultaneously at ±Δ (or at Δ = 0) eliminates mixing between Zeeman and Dipolar components of s2 (this is the ihMT effect); (iv) M s2 D relaxes to a thermal equilibrium value of zero.

| MRI signal simulations
The system of differential equations described above, referred to here as the Bloch-McConnell-Provotorov (BMP) equations, allow for MRI signal modeling of this system with arbitrary pulse sequences. Here, we consider rapid gradient echo sequences and make the common approximation that RF pulses act instantaneously. Under this approximation, the evolution during each pulse is Ṁ = ⟨ ⟩M where ⟨ ⟩ is the time average of over pulse duration , with solution M (t + ) = RM(t) where R = exp (⟨ ⟩ ). The matrix exponential preserves the block diagonal form of such that R consists of the normal rotation matrix for the free water in the upper left block, and a saturation matrix for the semisolid pool according to the mean square B 1 of the pulse in the lower right block; this is exactly the approach taken for modeling of "classic" MT in rapid sequences. 15 Evolution in the absence of RF is governed by Ṁ = M + C, which has the general solution M (t + TR) = SM (t) + (S − 1) −1 C for time period TR, where S = exp ( TR). It is hence possible to write the following signal expressions for SPGR: and bSSFP: The matrix accounts for spoiling or phase cycling in each case. For SPGR, we assume perfect spoiling of transverse magnetization and so = diag[0 0 1 1 1 1] is used to remove transverse components after each TR; for bSSFP, we assume phase cycling between successive pulses, hence = diag[ − 1 − 1 1 1 1 1]. The expression for SPGR gives the magnetization immediately after an RF pulse, while for bSSFP it is at the midpoint between 2 RF pulses. Equations 5 and 6 are generalizations of standard steadystate equations used for these sequences and give the same predictions as M s 0 → 0. Given that dipolar relaxation rates (R s 1D ) are often very fast (time constants <10ms), we have also examined the validity of treating RF pulses instantaneously by comparing with full time integration of the BMP equations.

| Time integration of BMP equations
Time integration for these sequences is cumbersome because hundreds of TR periods are required to reach a steady state, so a more direct method was devised. Equation 1 is first written as Ṁ = AM + C, and then reformulated as: where "1" is a scalar value, and "0" represents a row vector of zeroes matching the dimensionality of A and C ([0 0 0 0 0 0 0] in this case). The advantage of writing in this way is that the differential equation is homogeneous and the solution at any small time increment during which Ã is constant is M (t +Δt) = exp(Ã (t) Δt)M(t). This may now be integrated over one TR period, with a small time-step Δt allowing the shape of the RF pulse to be accounted for. The result is a single matrix that describes the evolution, such that: where n is an index and N = TR∕Δt. The product must be evaluated by left multiplication of successive matrices. Phase cycling or forced "perfect" transverse spoiling can be added by incorporating matrix into the product. The steady state solution requires that M is invariant after multiplication by X -i.e., it is the eigenvector of X with an associated eigenvalue of 1.
To remain consistent with Equations 5 and 6, the integration ̇M =ÃM Illustration of the multiband excitation pulses used in this work. All are spatially nonselective and have the same on-resonance contribution (shaded green on power spectra), which is responsible for flipping the free water magnetization. The off-resonant components illustrated have the same total power (red shaded area) and this can be split either between 2 symmetric bands, or 1 band with either a positive or negative frequency offset. Time domain profiles show that the 3-band pulses have the higher peak B 1 , as expected from Equation 10C is performed starting at the end of an RF pulse for SPGR, and starting at the middle of a TR period for bSSFP. As outlined above, semi is evaluated by considering the time domain representation of the individual frequency bands b{Δ j ;t} that compose the multiband pulses.

| METHODS
In this article we explore the use of nonselective multiband RF pulses (as proposed by Teixeira et al 10 ) with either 2 or 3 bands, to generate ihMT contrast in 3D gradient echo sequences. In all cases, the pulses have an on-resonance band responsible for flipping the free water magnetization; the 2-band pulses have a single off-resonance saturation band, whereas 3-band pulses have symmetric off-resonant bands (see Figure 3). Using these pulses, steady state solutions (Equations 5 and 6) are compared with time integration of the BMP equations (see Section 2.3), and expected contrast levels are explored for a range of tissue parameters taken from literature. Predictions are compared with phantom measurements, and in vivo images are also presented. Matlab code required for RF pulse generation and signal simulations has been posted at https ://github.com/mriph ysics/ ihMT_stead ystate, (hash 0db78e5 was used for the presented results). Phantom measurements are also available at the same site.

| Multiband RF pulse design
Multiband pulses with power spectra illustrated by Figure 3 can be generated by multiplying single-band on-resonance pulse b{0; t} by modulation function mod (t). A key parameter for quantifying saturation power is the mean square B 1 (⟨B 2 1 ⟩) 6,7 because this (or the root mean square, B rms 1 ) is a parameter easily accessible on many clinical MR systems. An SPGR or bSSFP sequence using single band excitation pulse b{0; t} has mean-square B 1 given by where is the pulse duration. This can be turned into the required multiband pulse using modulation functions defined as follows: where quantifies the ratio of off-resonant to on-resonant power, calculated by and ⟨B 2 1 ⟩ TOT is the desired total mean square B 1 such that the whole sequence B rms 1 ≡ √ ⟨B 2 1 ⟩ TOT . These expressions allow flexible design of multiband pulses that reach a required target total B rms 1 within a given pulse sequence; this depends on the properties of the pulse but also the sequence TR. All presented data used Gaussian pulse shapes for the single band pulses with time-bandwidth-product 2.26. For brevity, we will refer to single band pulses as 1B, 2 bands as 2B (or 2 + B/2 − B if distinction between sign of offset is important) and 3 bands as 3B. Function gen_MB_pulse.m in the accompanying code is an example implementation for generating these.

| Numerical simulations
Signal simulations were performed for white matter using parameters based on measurements of human internal capsule taken from Mchinda et al 7 relevant for 1.5T imaging: was taken from Sled and Pike. 9 A Super-Lorentzian absorption line-shape 14 was assumed for white matter. Because on-resonance absorption is included in the model, the on-resonance singularity in the Super-Lorentzian function was handled by extrapolating between ±1 kHz as proposed by Bieri and Scheffler. 16

| ihMT contrast from multiband steady-state imaging
Signal simulations were performed for short TR sequences using 1B, 2B, and 3B pulses and a range of other sequence parameters ( , TR, B rms 1 , flip angle, Δ). The effect of longterm dipolar order was assessed by modifying the white matter parameters to include a case with = 0 (no "dipolar" pool) and = 1 (whole semisolid pool contains "dipolar" component). The "ihMT difference" is defined as ; similarly the MT ratio is defined To account for asymmetry induced by chemical shift of the semisolid line, the definition is used for analysis of in vivo data. This is not necessary for numerical simulations because the line is centered on the water resonance in simulations. Because the single band images S 1B are used as a reference for calculating both MTR and ihMTR, it is important to note that these images will also contain on-resonance MT effects (e.g., see Gloor et al,) 17 although the effect is small (10b) 2 bands with offset −Δ: ) 3 bands with offsets ±Δ: for the sequences used in this work. This is inherently included in our calculations because the sum in Equation 4 includes the on-resonance band (i.e., Δ = 0). In general, the contrast in the S 1B image is influenced by MT and obviously also by free water relaxation effects. It should not necessarily be expected that MTR and ihMTR measured by steady-state sequences would be directly comparable to those obtained by other methods.

| Instantaneous RF approximation
The sequences used in this work typically used RF pulse durations that are not short with respect to the TR. Hence, the steady-state signal expressions Equations 5 and 6 were compared with time integration using the method described in Section 2.3. Numerical integration used step size Δt = 10 μs.
Simulations were performed for the white matter tissue parameters, except that T s 1D ≡ 1∕R s 1D was varied logarithmically from 0.1 ms to 10 ms. Balanced SSFP and SPGR sequences with TR = 5 ms and B rms 1 = 5 T were simulated for 1B, 2B, and 3B pulses for durations varying between 0.4 ms and 2.5 ms.
Bieri and Scheffler 18 proposed a correction to the "instantaneous approximation" bSSFP signal model for finite duration RF pulses, which amounts to a correction to R 2 to account for the fact that magnetization does not spend the full TR period in the transverse plane. This correction was also implemented and compared with full time integration. The original correction method is for a single pool model and amounts to modifying R 2 as a function of the true R 1 and R 2 and sequence properties. For our multi-compartment model, we use the same formulation but substitute R 1,2 for R f 1,2 in the expressions derived by Bieri and Scheffler. 18

| Phantom experiments
A phantom consisting of 4 sample tubes immersed in a bath of water for susceptibility matching was imaged on a Philips (Best, Netherlands) 1.5T Ingenia MR system. The samples included: water doped with MnCl 2 (~0.01 mM); chemically cross-linked bovine serum albumin (BSA, 10% w/w in water prepared as in Koenig et al 19 ); prolipid 161 (PL161, 15% w/w; Ashland Inc, Covington, KY; prepared as in Swanson et al 2 but excluding any T 1 reducing agent); and off the shelf hair conditioner (HC, TRESemmé, Unilever PLC, London, UK). The MnCl 2 phantom is not expected to exhibit any MT effect; BSA is a model substance for human tissue MT effects, 19 but is not expected to exhibit ihMT contrast because of protein cross-linking. PL161 2,3 and HC 8,20 have both been shown to exhibit strong MT and ihMT contrast as they contain lamellar liquid crystal structures.
The phantom was imaged using a series of bSSFP and SPGR sequences, all with TR = 5 ms, echo time = 2.5 ms, = 2.2 ms, Δ∕2 = 8 kHz, isotropic resolution 1.5 × 1.5 × 1.5 mm 3 , no parallel imaging acceleration, time per scan 57 s. A total of 64 images were obtained using flip angles 10° to 80° (steps of 10°) for bSSFP, and 2° to 16° (steps of 2°) for SPGR, and then for each of these using 1B, 2 + B, 2 − B, and 3B excitation pulses. All multiband pulses were computed for B rms 1 = 4.2 T, whereas the single band pulses used only the on-resonance component of these pulses, and so had variable B rms 1 over the range of flip angles. The data were analyzed by averaging the signal from each sample tube over a range of 45 mm along the length of each tube; segmentation was performed automatically using k-means clustering using all images as the feature space.
The signal models (Equations 5 and 6) were then fitted to these averaged measurements to determine best-fit model parameters for each sample. For each sample, data from all 64 images were fitted together as a single model fit with a common scaling factor added to account for unknown overall image scaling. Signals S 2 + B and S 2 − B were averaged before fitting because the model does not distinguish these. For SPGR sequences, Equation 5 was multiplied by an additional factor e −TE R f 2 to account for transverse relaxation at the echo time, here using the approximation that R For bSSFP sequences, the Bieri-Scheffler correction for finite pulse duration was applied. Matlab's fmincon optimization function was used with the interior point algorithm selected. Upper and lower bounds were used to constrain fits to plausible values; in each case, a range of bounds and initial starting points were tried out and the best fitting solutions in terms of normalized root mean square error (NRMSE) are reported. Uncertainties were computed by residual bootstrapping using 500 re-samplings for each fit. All phantom data were fit using a Gaussian lineshape function, because this was empirically found to yield the lowest NRMSE in all cases. A full set of fitting bounds used for generating the presented data is given in Supporting Information Table S1, which is available online.

| In vivo imaging
A single adult normal volunteer (age 24 years, male) was imaged on the same Philips 1.5T MR system, after giving written informed consent in line with local governance rules. Balanced SSFP images were acquired with flip angle 30°, TR = 5 ms, echo time = 2.5 ms, = 2.3 ms,Δ∕2 = 8 kHz, isotropic resolution 1.5 × 1.5 × 1.5 mm 3 , matrix 151 × 143 × 93, 9 signal averages, no parallel imaging acceleration, time per scan 9 min. A standard Philips 32 element head receiver was used. Separate images were acquired with 1B, 2 + B, 2 − B, and 3B excitation pulses. B rms 1 was 4.73 T for all multiband sequences and 0.76 T for single band; hence, the offresonance B rms 1 for multiband sequences was 4.67 T because the RMS values add in quadrature. Images were registered and brain extracted using FSL FLIRT and BET tools, respectively, 21,22 before MTR and ihMTR were calculated. MTR and ihMTR images were smoothed with 3D Gaussian kernel, standard deviation 1 mm. Figure 4 shows the predicted signals for the proposed sequences with TR = 5 ms, = 2 ms, Δ∕2 = 7 kHz, B rms 1 = 5 T over a range of flip angles (as described above the single band pulses have variable B rms 1 ). The "white matter" tissue parameters were used except that was set to 0 and 1 to examine extreme cases. For = 0, we see no difference between 2 or 3 band RF excitation; both lead to a strong reduction of signal compared with single band excitation, shaded blue. When = 1, signal S 3B is clearly lower than S 2B . This is expected because the 3B excitation decouples the dipolar ordered pool, leading to more efficient saturation of the semisolid Zeeman magnetization. Supporting Information Figure S1 plots ΔihMT and ihMTR as functions of flip angle for the = 1 case, showing that ihMTR is predicted to be in the order of 10% for bSSFP at lower flip angles, and similar for the SPGR around the Ernst angle. Figure 5 summarizes the comparison between the instantaneous pulse approximation and full integration methods described in Section 3.2.2. The left panels show example signal curves for 2 band excitation with B rms 1 = 5 T, TR = 5 ms, = 2 ms, Δ∕2 = 7 kHz, using white matter tissue parameters. Results are very similar for SPGR, but systematic differences exist for bSSFP that are resolved with application of the Bieri-Scheffler correction. The right-hand panels of Figure 5 explore the effect of making the instantaneous assumption by displaying the percentage deviation between instantaneous and full integration methods, relative to the full integration result. For 1B and 3B, large deviations only exist for bSSFP without the Bieri-Scheffler correction. For 2B excitation differences are present for both SPGR (up to 4.8%) and corrected-bSSFP (up to 2.4%) in the case of short T s 1D (<1 ms) and long . The Bieri-Scheffler correction was hence F I G U R E 4 Simulated signal curves for bSSFP and SPGR sequences for the model white matter tissue parameters, with set to either 0 or 1. When = 0 there is no dipolar ordered fraction and so signals from 2B and 3B pulses (S 2B & S 3B ) are expected to be the same. However, both are significantly attenuated in comparison with single band signal S 1B ; this is the classic MT effect, shaded blue on the plots. When = 1 the whole semisolid compartment is coupled to a dipolar fraction and there is a difference between S 2B and S 3B in this case of around 10% of S 1B ; this is the inhomogeneous MT effect, shaded green applied to all other bSSFP signal predictions presented in this study including Figure 4 Figure 6 shows the phantom imaging results with MTR and ihMTR calculated for single example flip angles from each image series (30° for bSSFP, 6° for SPGR). The BSA, HC, and PL161 samples all show an MTR of up to 50% for both sequences, while the water bath and MnCl 2 tube do not. The HC and PL161 samples additionally show an ihMTR of around 15%, suggesting they exhibit dipolar order effects. Figure 7 shows the model fits to data with the fitted parameters given in Table 1. The MnCl 2 phantom is expected to show no difference between any of the excitation pulse types and the best fit was obtained with M s 0 fixed at zero. Small discrepancies in the bSSFP data can be attributed to offresonance distortion from an air bubble in this sample, visible on Figure 6.

| Phantom experiments
The BSA phantom showed no large difference in signals from 2 + B, 2 − B, or 3B excitation, suggesting the absence of both ihMT effects and MT asymmetry; the best fit was achieved when constraining to be fixed at zero. For PL161, the best fit was obtained when was fixed to 1; however, for HC, a slightly improved fit was found when was allowed to vary. Uncertainty analysis also showed that some parameter estimates for HC are less precise than the others; this finding reflects instability of the model fitting, but also it should be noted that the residual bootstrapping method used is less valid for the situation with nonrandom structure in the residuals. Supporting Information Figure S2 plots MTR and ihMTR derived from these data.

| Choice of imaging parameters for in vivo experiment and imaging results
Results from phantom experiments indicated that similar ihMTR can be obtained using both bSSFP and SPGR, but contrast-to-noise ratio is higher for the former. We, therefore, focused on bSSFP imaging in vivo and sought to optimize acquisition parameters for the white matter tissue parameters listed in Section 3.2. The steady-state method has some particular properties: (i) the optimal flip angle to use depends on the steady-state of the whole system, not just the dipolar component; (ii) hardware limits constrain the achievable parameters but sometimes counterintuitively. Figure 8 illustrates some of these trade-offs. For example, Figure 8A shows that the peak ΔihMT occurs at flip angle ~30° and offset ~8 kHz for = 2 ms and B rms 1 = 5 T, which is the maximum allowed within SAR constraints on the MR system used here. Figure 8B shows that increasing the B rms 1 would continue to increase the contrast. Parts C and D show ihMTR instead; here the peak value occurs at a lower flip angle but this is because the reference single band image has a lower signal at that flip angle. Setting the parameters to maximize ΔihMT maximizes the effect size relative to the total magnetization; ihMTR in the region of 5% is then expected for white matter. Figure 8E shows that use of single multiband pulses for both saturation and excitation leads to a coupling between hardware constraints; to achieve a certain B rms 1 for a fixed the peak B 1 needed (i.e., the peak amplitude of the RF pulse) increases as TR increases. This is because reducing TR leads to an increase in B rms 1 for a given RF pulse, as the same energy F I G U R E 7 Fits to data (details in text, fitted parameter values in Table 1). In each case, all acquired datapoints were fitted to the model with a single parameter set. Excellent fits were obtained for water, BSA, and PL-161 but for HC some error was seen in the SPGR fit. is delivered at a higher rate. Increasing pulse duration reduces this maximum, but also reduces B rms 1 and may necessitate an increase in readout bandwidth, affecting signal-to-noise ratio. The human imaging protocol was designed based upon these considerations. Figure 9 shows the acquired human brain MTR and ihMTR images. As predicted the ihMTR in white matter is close to 4-5%, while MTR is nearer to 30%. The images show visibly different contrasts with ihMTR appearing more specific to white matter, and with the corticospinal tracts clearly visible.

| DISCUSSION
This work has demonstrated inhomogeneous MT contrast in steady-state gradient echo sequences using multiband pulses. Typical ihMT sequences use a saturation preparation step followed by a readout step. The combined multiband approach instead uses the same RF pulse to both excite free water and provide saturation of semisolid magnetization. This study presents theory for modeling the resulting signals, designing the RF pulses, and experimental validation in phantoms and in vivo. The behavior of bSSFP and SPGR sequences is exemplified by Figure 4. As expected, 3B pulses (on-resonance and dual off-resonance) lead to a lower signal than 2B (on-resonance and single off-resonance) when a dipolar order term is present despite containing the same overall power. This is because the 2B pulses populate the off-diagonal entries to semi in Equation 4 while 3B pulses close this pathway, leading to more efficient saturation of the Zeeman semisolid magnetization. Human brain imaging results ( Figure 9) demonstrate similar contrast to other published ihMTR measurements (e.g., see references 4,6,7,23).   The derived signal equations (Equations 5 and 6) predict ihMTR in the region of 4-5% in white matter using literature tissue parameters for our sequences, and this correlates well with observations. The mathematical signal model used in this work was developed by Varma and co-authors over many publications. 1,8 In addition to explicitly including a sum over frequency components necessary for handling multiband pulses, the contribution of this study is to assemble these into a matrix formalism that also includes transverse components of the free water pool, and to then demonstrate that the same matrix equations used to model steady-state properties in simpler systems can also be used to model ihMT.
In particular, the commonly used instantaneous pulse approximation was used to simplify the calculations, and accuracy was investigated by comparing with full time integration. It may be expected that this approximation would fail for ihMT because the dipolar relaxation times are typically of the order of a few milliseconds, not much greater than the RF pulse durations or indeed the TRs examined. However, we found that the approximation does hold for SPGR sequences in general, with appreciable differences only when 2B pulses are used with duration ≳ T s

1D
. The fact that differences only occur with 2B pulses is a clear indication that this is due to the dipolar terms, because these are excluded for 1B and 3B. It is also reasonable that errors would occur when is not small compared with T s 1D because we would then not properly account for relaxation of the dipolar pool during saturation. For bSSFP, there are much larger discrepancies associated with the instantaneous approximation, however, these largely vanish when the Bieri-Scheffler correction 18 is applied to R f 2 to account for the fact that the magnetization spends less than the full TR period in the transverse plane. Although derived for a standard Bloch equation model, it is perhaps interesting that this correction still holds here. An explanation for this may come from the original study, 18 which showed that the longitudinal term is only marginally affected. The additional terms in our MT/ihMT model can be thought of as "longitudinal," because they couple only to the free water longitudinal magnetization and have no transverse components; hence, it is reasonable that they would behave in the same way.
Another contribution of this work is the development of the efficient eigenvector based direct integration method. The method is similar to other approaches proposed for quantitative MT 24 or CEST, 25 but to our knowledge is unique in directly computing the steady-state signal for rapid imaging sequences in this way and could be generalized to other multi-pool (or single pool) models. Our results indicate that the instantaneous approximation does hold for the relatively short RF pulses used in this work, but it is not guaranteed to for longer saturation pulses. Portnoy and Stanisz 26 proposed a method to address this issue, which works by integrating Equation 1 through time for specific pulses. A problem is F I G U R E 9 In vivo image results showing MTR (A) and ihMTR (B) that because this equation is not homogenous, the effect of a long RF pulse cannot be easily written as a compact product of matrices, as there are multiple terms to consider. The "augmented" matrix approach (i.e., Equation 7) provides a way to sidestep this issue to make the effect of any long pulse simply the product of many matrices. Eigenvector decomposition can then be used to obtain the steady-state directly without forward integration over multiple TR periods. The relative acceleration afforded by direct calculation is proportional to the number of TR periods that would otherwise have to be simulated; this is typically in the hundreds. The high efficiency of the eigenvector based approach could make it a viable simulation route in cases that the instantaneous approximation does not hold.
Phantom experiments showed that this model can achieve good fits to data, with PL161 the most relevant model for ihMT also studied by others 2,3 an excellent fit to data was achieved with NRMSE = 0.92%. The fitted T s 1D = 20.7 ms agrees with that of Swanson et al, 2 other quantities are not available in literature or are not comparable as our sample was prepared without a T 1 shortening agent; however, all values appear plausible. Of interest, the estimated T s 1Z was in the 100-to 200-ms range for the phantoms with MT effects, which is quite far from the typically assumed 1 s used by many in the literature. Figure 7 also shows that there is almost no difference between 2 + B and 2 − B for PL161, suggesting very little asymmetry of the absorption line. The fits for HC have slightly larger NRMSE (2.4%), primarily due to a poor fit for the SPGR measurements; a potential explanation for this is that this substance may contain extra unmodeled magnetization pools. We generally observed a high degree of degeneracy in fitting for all phantoms with MT effects. In particular, this was true for the T 1 of free and semisolid compartments, exchange rate, and fraction . We hence limit our conclusions only to the fact that the model does fit the observed data well with plausible parameters, and do not claim that this experiment constitutes a proper quantitative measurement of these parameters. An obvious extension to better explore the lineshape properties would be to add measurements at different frequency offsets, for example. Note that, while some parameters were fixed for fitting, the reported solutions were the ones that led to the lowest NRMSE in each case (i.e., as far as could be determined, allowing fixed parameters to vary did not lead to reduced NRMSE).
The use of multiband pulses in this study was motivated by related work that showed that this type of RF excitation can stabilize gradient echo based relaxometry by equalizing MT effects between measurements. 10 Teixeira et al have shown that estimated T 1 is a function of B rms 1 (see fig. 7 in Teixeira et al 10 ). An extension to that work incorporating the theory developed here will be to investigate whether changes in apparent T 1 when using 2B and 3B pulses can be related to ihMT parameters. Similar "pseudo-quantitative" work from Geeraert et al 27 suggests this is a possibility.
In addition to the potential for steady-state relaxometry, rapid gradient echo sequences, particularly bSSFP, offer high signal-to-noise ratio per unit time and can generate high B rms 1 using relatively short RF pulses because these are repeated rapidly with short TR. The ihMT effect is small; Figure 8 suggests that the contrast in absolute terms should be approximately 0.7% of the total available magnetization (M 0 ), for example. This is limited by B rms 1 , which increases directly with SAR; indeed, one motivation for conducting experiments at 1.5T is that these scanners can typically access much larger B rms 1 than higher field systems. Measured ihMTRs by our method are of the order of a few percent, which are lower than those measured by other methods, 7 although this might largely be explained by a factor of 2 difference between the definition of ihMTR used in this study and in other works (see for example, Girard et al 6 ). Regardless, the measurements are not directly comparable to those from other methods because ihMTR depends also on the signal value in the 1B image, which depends on the relaxation properties of the free water pool, and to a smaller extent on-resonance MT effects. A direct comparison of efficiency between steady-state and other existing approaches will constitute future work. Proposed methods to boost contrast-tonoise ratio, such as cycling to higher RF power when the low k-space frequencies are sampled 7 or the observation that low duty-cycle pulsed saturation 23 can enhance contrast, will also be explored for the multiband gradient echo sequences used in this work.

| CONCLUSIONS
Steady-state imaging sequences generating inhomogeneous MT contrast have been demonstrated using spatially nonselective multiband RF excitation pulses. A matrix-based approach for quantifying steady-state signals was developed and demonstrated to provide good fits to phantom data, and to make good predictions for observed contrast in in vivo human brain images using balanced SSFP.