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 myelin4 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 single‐frequency saturation pulses have been proposed to generate ihMT weighted images6, 7 demonstrating strong white‐matter specificity. More recently, quantitative approaches to measure the dipolar relaxation time8 and influence of ihMT on free water T1 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 sequences9 often consist of separate saturation and excitation pulses preceding a gradient echo readout (Figure 1A), while a recently proposed ihMT sequence7 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.

Figure 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 al7 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

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 pulses are used. We investigate the validity of this matrix‐based model, compare predictions with phantom experiments, and present in vivo ihMT ratio images formed using the proposed sequences.

THEORY

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 (Zeeman) magnetization of the free water pool interacts with (Zeeman) magnetization of the semisolid pool, written . The semisolid pool has a very short transverse relaxation time (), so transverse components are typically excluded. The ihMT effect arises in some semisolids because of the presence of a dipolar interaction term (D) as well as Zeeman term (Z) in the quantum Hamiltonian. According to Provotorov theory as presented by Goldman,11 this may be described by subdividing the net magnetization in the semisolid pool between Zeeman and dipolar ordered pools and , respectively; here we use the formalism adopted by Lee et al12 to represent both terms as dimensionless polarizations. The magnetization of this system can be written as . The model has been further updated by Varma et al1 to include multiple semisolid pools, some of which do not contain a dipolar order term. Varma et al1 showed that a model with two semisolid pools, one without and one with dipolar order (denoted s1 and s2, respectively) provides a good fit to experimental data. The overall state is, therefore, described by vector whose temporal evolution is governed by:

Matrices and describe evolution in the absence of RF:where are the free water relaxation rates (inverse of relaxation times); and 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 and 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 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 and , respectively, of the total semisolid compartment and we adopt the convention . Figure 2 illustrates the relationships between compartments in this model.

Figure 2

Multi‐compartment model consists of free water in contact with a semisolid compartment that has two sub‐parts 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 is the (angular) frequency offset of the j th band and is the time domain pulse waveform associated with component frequency . It then follows that matrix describing interactions due to RF pulses may be written as:

is simply the Bloch equation for the free pool. For the multiband pulses used in this work, only the on‐resonance component need be considered. describes absorption by the semisolid pools and is explicitly written as a sum over frequency bands as suggested by Goldman11 for the case of “fast modulation”. The RF saturation rate is defined as , where is the absorption lineshape of the semisolid and is its transverse relaxation time. The term 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 . The term 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 is exactly as commonly used to describe “classic” MT effects; (ii) exchange between and is purely RF driven, while exchange between and occurs continuously; (iii) coupling terms in are linear in so equal off‐resonance RF power applied simultaneously at (or at ) eliminates mixing between Zeeman and Dipolar components of s2 (this is the ihMT effect); (iv) 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 where is the time average of over pulse duration , with solution where . The matrix exponential preserves the block diagonal form of such that 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 B1 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 , which has the general solution for time period TR, where . 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 is used to remove transverse components after each TR; for bSSFP, we assume phase cycling between successive pulses, hence . 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 steady‐state equations used for these sequences and give the same predictions as . Given that dipolar relaxation rates () 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 , 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 . This may now be integrated over one TR period, with a small time‐step 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 . 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 is invariant after multiplication by — i.e., it is the eigenvector of with an associated eigenvalue of 1. To remain consistent with Equations 5 and 6, the integration 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, is evaluated by considering the time domain representation of the individual frequency bands that compose the multiband pulses.

METHODS

In this article we explore the use of nonselective multiband RF pulses (as proposed by Teixeira et al10) 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/mriphysics/ihMT_steadystate, (hash 0db78e5 was used for the presented results). Phantom measurements are also available at the same site.

Figure 3

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 B1, as expected from Equation 10C

Multiband RF pulse design

Multiband pulses with power spectra illustrated by Figure 3 can be generated by multiplying single‐band on‐resonance pulse by modulation function . A key parameter for quantifying saturation power is the mean square B1 ()6, 7 because this (or the root mean square, ) is a parameter easily accessible on many clinical MR systems. An SPGR or bSSFP sequence using single band excitation pulse has mean‐square B1 given bywhere 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 byand is the desired total mean square B1 such that the whole sequence . These expressions allow flexible design of multiband pulses that reach a required target total 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 al7 relevant for 1.5T imaging: , , , , , , ; was taken from Sled and Pike.9 A Super‐Lorentzian absorption line‐shape14 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, , flip angle, ). The effect of long‐term dipolar order was assessed by modifying the white matter parameters to include a case with (no “dipolar” pool) and (whole semisolid pool contains “dipolar” component). The “ihMT difference” is defined as where is the signal from the 2B acquisition, etc. The “ihMT ratio” is then defined as ; similarly the MT ratio is defined as . 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 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 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., ). In general, the contrast in the 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 . Simulations were performed for the white matter tissue parameters, except that was varied logarithmically from 0.1 ms to 10 ms. Balanced SSFP and SPGR sequences with TR = 5 ms and were simulated for 1B, 2B, and 3B pulses for durations varying between 0.4 ms and 2.5 ms.

Bieri and Scheffler18 proposed a correction to the “instantaneous approximation” bSSFP signal model for finite duration RF pulses, which amounts to a correction to R2 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 R2 as a function of the true R1 and R2 and sequence properties. For our multi‐compartment model, we use the same formulation but substitute for 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 MnCl2 (~0.01 mM); chemically cross‐linked bovine serum albumin (BSA, 10% w/w in water prepared as in Koenig et al19); prolipid 161 (PL161, 15% w/w; Ashland Inc, Covington, KY; prepared as in Swanson et al2 but excluding any T1 reducing agent); and off the shelf hair conditioner (HC, TRESemmé, Unilever PLC, London, UK). The MnCl2 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. PL1612, 3 and HC8, 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, = 8 kHz, isotropic resolution 1.5 × 1.5 × 1.5 mm3, 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 , whereas the single band pulses used only the on‐resonance component of these pulses, and so had variable 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 and were averaged before fitting because the model does not distinguish these. For SPGR sequences, Equation 5 was multiplied by an additional factor to account for transverse relaxation at the echo time, here using the approximation that . 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, = 8 kHz, isotropic resolution 1.5 × 1.5 × 1.5 mm3, 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. was for all multiband sequences and for single band; hence, the off‐resonance for multiband sequences was 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.

RESULTS

Expected signal curves and instantaneous pulse approximation

Figure 4 shows the predicted signals for the proposed sequences with TR = 5 ms, = 2 ms, = 7 kHz, over a range of flip angles (as described above the single band pulses have variable ). 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 is clearly lower than . 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 and ihMTR as functions of flip angle for the 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 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 (S2B & S3B) are expected to be the same. However, both are significantly attenuated in comparison with single band signal S1B; 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 S2B and S3B in this case of around 10% of S1B; this is the inhomogeneous MT effect, shaded green

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 , TR = 5 ms, = 2 ms, = 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 (<1 ms) and long . The Bieri‐Scheffler correction was hence applied to all other bSSFP signal predictions presented in this study including Figure 4.

Figure 5

Examination of the instantaneous pulse approximation compared with full time integration. Left, Example SPGR and bSSFP curves for , TR = 5 ms, = 2 ms, = 7 kHz for 2B pulses. The bSSFP curve without Bieri‐Scheffler correction diverges from the full integration curve at higher flip angles. Right, Percentage deviation between instant approximation and full integration for 1B, 2B, and 3B pulses with different and , flip angle = 60° for bSSFP and 7.1° (Ernst angle) for SPGR. Larger deviations are observed for 2 band pulses with short and longer , and bSSFP without Bieri‐Scheffler correction

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 MnCl2 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 MnCl2 phantom is expected to show no difference between any of the excitation pulse types and the best fit was obtained with fixed at zero. Small discrepancies in the bSSFP data can be attributed to off‐resonance distortion from an air bubble in this sample, visible on Figure 6.

Figure 6

Phantom experiment results. Left, bSSFP image for 30° flip angle (1B) showing the phantom. Green arrow is marking the position of an air bubble in the MnCl2 phantom. Right, images of MTR and ihMTR as defined in the text, for 30° flip angle for bSSFP and 6° flip angle for SPGR. All MTR and ihMTRs are quoted in percentage units using the relevant single band image as a reference. Notably, the BSA phantom has a strong MTR but no ihMTR, while the PL161 and HC show both effects

Figure 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. BSA shows no ihMT effect (2B & 3B identical), and in all cases 2+B and 2−B are very similar, suggesting little to no MT asymmetry in any sample. Signal values are as reconstructed by the scanner (arbitrary units) but are consistent between plots

Table 1

Fitted parameters for phantom experimentsa

	T1f (ms)	T2f (ms)	M0s	T1Zs (ms)	T1Ds (ms)	T2s (µs)	k (s‐1)	δ
MnCl2	2174 ± 152	683 ± 45	0 (fixed)	–	–	–	–	–
BSA	1643 ± 29	53.2 ± 0.1	0.075 ± 0.001	130 ± 0	–	17.6 ± 0.0	70.1 ± 0.0	0 (fixed)
PL161	1912 ± 80	72.8 ± 0.4	0.150 ± 0.003	202 ± 2	20.7 ± 0.0	17.5 ± 0.1	46.7 ± 0.0	1 (fixed)
HC	1975 ± 127	158 ± 4	0.060 ± 0.004	189 ± 6	23.3 ± 0.1	19.2 ± 0.1	67.8 ± 0.1	0.683 ± 0.088

Note that the relaxation rates are reported here as times (inverse of rate) because these are more intuitive; they were, however, fitted as rates. Uncertainties were estimated using residual bootstrapping.

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 occurs at flip angle ~30° and offset ~8 kHz for = 2 ms and , which is the maximum allowed within SAR constraints on the MR system used here. Figure 8B shows that increasing the 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 maximizes the effect size relative to the total magnetization; ihMTR in the region of 5% is then expected for white matter.

Figure 8

Results of contrast simulations for white matter. A, as a function of and flip angle, for fixed TR = 5 ms, = 2 ms, and . B, as a function of and for fixed TR = 5 ms, = 2 ms, and flip angle 30°. C,D, ihMTR computed from (A) and (B), respectively. E, Maximum B1 (i.e., peak B1 in each RF pulse) as a function of TR and for fixed = 2 ms and flip angle 30°. The shaded regions show inaccessible solutions on our 1.5T system. RF duty cycle limits suggest and SAR limits mean . In addition, max B1 cannot exceed 20

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 for a fixed the peak B1 needed (i.e., the peak amplitude of the RF pulse) increases as TR increases. This is because reducing TR leads to an increase in for a given RF pulse, as the same energy is delivered at a higher rate. Increasing pulse duration reduces this maximum, but also reduces 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.

Figure 9

In vivo image results showing MTR (A) and ihMTR (B)

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 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 . 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 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 correction18 is applied to 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 MT24 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 Stanisz26 proposed a method to address this issue, which works by integrating Equation 1 through time for specific pulses. A problem is 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 others2, 3 an excellent fit to data was achieved with NRMSE = 0.92%. The fitted = 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 T1 shortening agent; however, all values appear plausible. Of interest, the estimated 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 T1 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 T1 is a function of (see fig. 7 in Teixeira et al10). An extension to that work incorporating the theory developed here will be to investigate whether changes in apparent T1 when using 2B and 3B pulses can be related to ihMT parameters. Similar “pseudo‐quantitative” work from Geeraert et al27 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 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 (M0), for example. This is limited by , which increases directly with SAR; indeed, one motivation for conducting experiments at 1.5T is that these scanners can typically access much larger 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 al6). 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‐to‐noise ratio, such as cycling to higher RF power when the low k‐space frequencies are sampled7 or the observation that low duty‐cycle pulsed saturation23 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.

FIGURE S1 Top for bSSFP and SPGR as a function of flip angle for the case plotted in Figure 4. Here is shown as a percentage of the overall M0 which was arbitrarily fixed as 1.0 in this work. Middle, ihMTR for the same data, here shown as a percentage of the signal from the equivalent single band image (i.e. % of S1B). Close to 10% is achieved for both bSSFP and SPGR at lower flip angles. Bottom, the multiband pulses used in this work were computed to maintain the total ; hence, increasing the on‐resonance flip angle results in a reduction of the total off‐resonance saturation power (i.e., ) that creates ihMT contrast. This plot shows for the flip angles used to make this figure (and Figure 4) showing that the effect is small for this sequence. Note that the pulse properties in this simulation were similar, although not exactly the same as those used in vivo. The reader should be aware that the on‐resonance (i.e. ) is larger than the gap between the curve and the line because the RMS values add in quadrature

FIGURE S2 MTR and ihMTR as a function of flip angle, derived from the data plotted in Figure 7; please note that these plots use separate y‐axes for MTR and ihMTR because of their different ranges. For the two ihMT phantoms (PL‐161 and HC), the MTR is much larger than ihMTR. BSA shows large MTR but small ihMTR; a small degree of ihMTR is seen for BSA at high flip angles with the SPGR data but this could not be explained easily by model fitting. Nonzero ihMTR observed in water was attributed to errors from off‐resonance related to an air bubble

TABLE S1 Upper and lower fit bounds (UB/LB) for each parameter

Click here for additional data file.