PRE-PRINT: SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE Magnetization Transfer-Mediated MR Fingerprinting

Purpose: Magnetization transfer (MT) and inhomogeneous MT (ihMT) contrasts are used in MRI to provide information about macromolecular tissue content. In particular, MT is sensitive to macromolecules and ihMT appears to be specific to myelinated tissue. This study proposes a technique to characterize MT and ihMT properties from a single acquisition, producing both semiquantitative contrast ratios, and quantitative parameter maps. Theory and Methods: Building upon previous work that uses multiband radiofrequency (RF) pulses to efficiently generate ihMT contrast, we propose a cyclic-steady-state approach that cycles between multiband and single-band pulses to boost the achieved contrast. Resultant time-variable signals are reminiscent of a magnetic resonance fingerprinting (MRF) acquisition, except that the signal fluctuations are entirely mediated by magnetization transfer effects. A dictionary-based low-rank inversion method is used to reconstruct the resulting images and to produce both semiquantitative MT ratio (MTR) and ihMT ratio (ihMTR) maps, as well as quantitative parameter estimates corresponding to an ihMT tissue model. Results: Phantom and in vivo brain data acquired at 1.5T demonstrate the expected contrast trends, with ihMTR maps showing contrast more specific to white matter (WM), as has been reported by others. Quantitative estimation of semisolid fraction and dipolar T1 was also possible and yielded measurements consistent with literature values in the brain. Conclusions: By cycling between multiband and single-band pulses, an entirely magnetization transfer mediated 'fingerprinting' method was demonstrated. This proof-of-concept approach can be used to generate semiquantitative maps and quantitatively estimate some macromolecular specific tissue parameters. PRE-PRINT: SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE


Magnetization transfer (MT) and inhomogeneous MT (ihMT) contrasts are used in MRI to
provide information about macromolecular tissue content. In particular, MT is sensitive to macromolecules and ihMT appears to be specific to substances with non-zero dipolar magnetization order, including the lipid bilayers that form myelin (1)(2)(3)(4). Recently, ihMT measurements have: compared favorably against other myelin imaging metrics (5,6); correlated strongly with fluorescence microscopy findings (7); and provided sensitivity for the assessment of demyelinating conditions, such as multiple sclerosis (8). MT and ihMT contrasts are usually obtained by sequences with off-resonance radiofrequency (RF) saturation pulses (that only affect the semisolid pool), followed by readout periods for measurement. For example, ihMT gradient echo methods have been developed and optimized to achieve high-resolution whole-brain imaging at 1.5T (9).
The ihMT effect arises because semisolid magnetization can be modelled as containing pools of both Zeeman and dipolar order that can be made to exchange by the presence of off-resonant RF irradiation; dual frequency saturation with equal and opposite frequency offsets cancels this interaction. In our previous work, it was shown that ihMT contrast can also be generated using non-selective multiband pulses that perform off-resonance saturation and on-resonance excitation simultaneously (10). These pulses were originally proposed to control for MT effects in variable flip angle techniques by ensuring constant RF power across all flip angles, resulting in more stable relaxometry measurements (11). In order to generate ihMT contrast, images using 2-band excitation (one on-resonance and one off-resonance band) leading to Zeeman-dipolar coupling, were compared with 3-band excitation (one on-resonance and two equal and opposite off-resonance bands) in which this coupling is cancelled. Use of the same total power results in the same classical MT effect, for tissues where dipolar order can be neglected.
Resultant ihMT ratios (ihMTRs) from this type of sequence are generally small, for example ~4% in white matter (WM), which is similar to other steady-state measurement techniques (9). In order to boost this contrast, Varma et al. recently proposed use of a low duty-cycle RF saturation scheme where saturation pulses are concentrated into short time periods with interleaved recovery time during which data are acquired (12). With this as motivation, in this work we propose a modulated cyclic-steady-state sequence employing a balanced steady-state free precession (bSSFP) acquisition with multiband pulses that are alternated over time to create a similar low duty-cycle saturation effect. Each RF pulse contains an on-resonance component, meaning that data is continuously read out, improving the efficiency of encoding.

PRE-PRINT: SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE
By holding the on-resonance flip angle constant while periodically alternating off-resonant saturation, the proposed sequence results in signal fluctuations that are entirely mediated by MT. This time-varying signal is then used for quantitative parameter estimation by applying current MR Fingerprinting (13) (MRF) methods. Images are spatially encoded using an undersampled tiny golden angle radial k-space trajectory (14) and reconstructed using dictionary-based low-rank inversion (LRI) (15) previously used for MRF. The proposed sequence is essentially an MRF acquisition but one which should yield only constant signals if tissues follow the standard Bloch Equations.
We refer to our approach as magnetization transfer-mediated MR fingerprinting (MT-MRF) to reflect these similarities and highlight the fact that induced signal variations are mediated by MT effects. As well as generating quantitative parameter estimates, we also show that the reconstructed time-series data can be used to produce semiquantitative MTR and ihMTR maps. We present an evaluation of the proposed method using phantoms, as well as MT-MRF brain images acquired from five healthy subjects. A comparison with a previous steady-state method is also included (10).

Sequence Description
The proposed sequence consists of a 3D bSSFP acquisition with constant flip angle (FA) and repetition time (TR), in which the excitation pulses are periodically switched between single-band and multiband pulses with different variants as shown in Figure 1. One cycle comprises a number, MB of 2-band pulses, followed by a block of 1B single-band pulses, a block of MB 3-band pulses, then another block of 1B single-band pulses (to allow semisolid saturation to reduce before the next cycle). The on-resonance bands of each pulse are identical, hence a tissue with no semisolid compartment would give constant signal. Tissues with a semisolid component respond transiently as they saturate strongly during multiband periods due to off-resonant bands and recover during single-band periods. Tissues with significant dipolar order effects will also respond differently to the 2-band and 3-band pulses. Example signal curves are shown in Figure 1. Figure 1: The proposed acquisition uses a rapid gradient echo sequence alternating between the pulses shown in the top row in blocks (bSSFP is used in this work); MB indicates the number of multiband pulses in a block and 1B is the equivalent number of single-band pulses. The on-resonance part of each pulse (with power α) is identical, so free water (no MT effect) would give constant signal throughout the acquisition. The switched off-resonance bands will affect semisolid proton saturation, resulting in modulation of the free water signal via MT as shown. 2-band (2B) and 3-band (3B) pulses are power-matched (with a combined off-resonance power β) and alternated to give ihMT contrast.

PRE-PRINT: SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE
The sequence repeats many times; the total number of excitations in one cycle (N cycle = 2 MB + 2 1B ) is typically of the order of 1000, such that the cycle repeats after a few seconds for the TRs (~5ms) achieved using bSSFP with multiband pulses. After a small number of repeats, the signal reaches a cyclic-steady-state whose MT or ihMT contrast is dictated by parameters such as 1B and MB as well as FA, TR, off-resonance power and offset frequency. Qualitative comparison of the efficiency of ihMT contrast generation for different parameter combinations is achieved by defining η as peak ihMT contrast generated per square root of time, where TR is the sequence repetition time. Times 1 and 3 are at the ends of 2B and 3B pulse periods, as in Figure 1.

Spatial Encoding and Image Reconstruction
The cycling of the acquisition sequence can repeat indefinitely and can be made independent of spatial encoding. Since the pulses are spatially non-selective, the encoding must be 3D and cover the entire field-of-view (FOV; the head). Full 3D characterization of each of the N cycle different contrast states is infeasible using regular encoding methods. Instead we take advantage of the fact that changes in signal are compressible in the temporal direction and can be described by a low-rank representation determined from a dictionary of simulations computed for a range of expected tissue properties. As shown by McGivney et al. (16), the dictionary can be compressed by singular value decomposition (SVD) and its low-rank approximation can be obtained by multiplication of the dictionary and the matrix of left singular vectors , truncated to a certain rank R to yield R ϵ ℂ TNC×RNC where T is the number of acquired timepoints, N is number of voxels and C is number of coils. The reconstruction problem was previously formulated for MRF as (15)(16)(17): Here, ϵ ℂ TKC is the k-space signal where K is k-space trajectory length. In Equation 2, ϵ ℂ RNC×RNC performs nonuniform fast Fourier transform and gridding operations; ϵ ℂ RNC×RN represents coil sensitivity maps; is a time-series of images sought to be reconstructed and ̃ is its low-rank approximation ( = R̃ so ̃= R , where ϵ ℂ TN and ̃ ϵ ℂ RN ). Patch-based regularization was added to the problem, as in Bustin et al. (18). constructs a 3D local tensor around voxel b by concatenating local (in each patch), non-local (between similar patches in a neighborhood) and contrast voxels along each dimension (19,20). This problem can be solved using the alternating direction method of multipliers (21); further details can be found in Ref. (18).
There is some flexibility in how much data needs to be collected and how this should be spatially encoded. It has been found that radial or spiral k-space encoding leads to better conditioning of the reconstruction problem than Cartesian sampling (22)(23)(24). Furthermore, the framework allows for a flexible amount of undersampling to be used. In total, the amount of data collected would be sufficient to reconstruct N V fully sampled volumes (ignoring temporal modulation of the signal), while the number of actually reconstructed volumes is R.

Signal Model and Efficient Simulations
The tissue model used in this work (1) consists of a single pool of 'free water' protons, denoted f, and a semisolid pool that contains two sub-compartments: one without dipolar order effects, denoted s1 (fractional size 1 − ), and the other with dipolar order, denoted s2 (fractional size ).

Magnetization is given by vector
T -subscript Z represents Zeeman ordered longitudinal magnetization and D denotes dipolar ordered longitudinal magnetization.
Malik et al. (10) proposed an efficient means for simulating the steady-state behavior of a sequence of RF pulses and delay periods. The governing equation is written in an 'augmented' form ̃̇=̃̃, whose solution for a period in which ̃ is constant is: ̃( + Δ ) = exp(̃( )Δ )̃( ).
Evolution over multiple periods can be summarized by a single matrix, which is the product of many matrix exponential terms. The resultant signal can be calculated by finding the eigenvector of this matrix whose eigenvalue is one. Specific details for MT-MRF are given in Supporting Information.

Sequence Design
Sequences of the type shown in Figure 1 were simulated using internal capsule tissue parameters from Mchinda et al. (9) and η values calculated for different parameter ranges, as shown in Figure 2. Several constraints were considered according to scanner hardware limitations: maximum pulse amplitude, B 1,max < 20μT, TR > 2τ (due to RF amplifier duty-cycle limits) and a 'timing' constraint, TR > 2.5 + τ where 2.5ms is readout duration and τ is RF pulse duration. Simulated signals were compared to our previous steady-state ihMT ('ss-ihMT') method with scan parameters matching the in vivo scheme used in Malik et al. (10). The "maximum contrast" solution in Figure 2 was found using a genetic algorithm (ga in MATLAB 2019a). The "scan" case was used for all experiments. 2-band pulse period, then 1-band etc.). ihMT contrast is the difference in signals at the end of the 2B and 3B periods.
'Scan' parameters generate greater contrast than the steady-state sequence. The 'maximum contrast' sequence obtains slightly higher contrast, while the 'forbidden' scheme (green) yields far higher contrast but violates peak B1 constraints.
The sequence design for MT-MRF was chosen to give large ihMT contrast, with the constraint that all on-resonance flip angles be constant, i.e. that all signal fluctuations are MT-mediated only.
To understand how well this sequence could be used to estimate underlying model parameters, we computed the Cramér-Rao Lower Bound (CRLB) (26)(27)(28)(29) and used this to estimate parameter-to-noise ratios (PNRs). Since the PNR and image signal-to-noise ratio (SNR) are both dependent on the underlying acquisition noise level, we have quoted the ratio PNR/SNR that independently quantifies the precision relative to image SNR. Figure   feature the most fixed parameters and the right-most columns feature the most estimated parameters). To achieve a reasonable PNR, only three parameters can be estimated simultaneously. Combination 2 is used for all dictionary fits.

Dictionary Generation, Low-Rank Basis and Parameter Estimation
Two types of dictionary were required for this work. Firstly, a 'full' dictionary denoted included variation of all model parameters. Having nine free parameters and a long temporal duration of each atom (1200 readouts) meant that the dictionary would quickly become too large for the memory of the PC used for calculation (8(16) × Intel® Core™ i7-5960X 3.00 GHz CPU, 64 GB RAM). Hence a version of was constructed using coarse sampling with the sole purpose being to discover a lower rank basis as detailed above. Secondly for parameter estimation, a reduced dictionary ̃ was constructed using finer sampling of the three estimated parameters (with all others fixed).
was created with ∽1 million atoms using parameter value ranges: 1 = 0. The low-rank basis R was constructed by performing an SVD on . A further reduction by randomly sampling 600,000 atoms was needed to perform the SVD; the random sampling was done ensuring an equal number of "no semisolid" ( = 0), "no dipolar order" ( ≠ 0, = 0) and "ihMT" entries ( ≠ 0, ≠ 0). The ability of the resulting basis to reproduce arbitrary tissue signals was assessed by projecting some test signals into the low-rank space via multiplication by R (truncated to different ranks) then back into the time domain by R . These approximated signals were compared to their ground-truth simulated equivalents ( Figure 4). Test signals were not themselves dictionary atoms but represented expected tissue values for cerebrospinal fluid (CSF), GM and WM. CSF was assigned: 1 = 3s, 2 = 2s, and = 0, whilst GM/WM took their respective values from Section 3.2.
Reduced dictionary ̃ was generated to estimate only 1 , f and 1D simultaneously using 100 increments of each parameter over their feasible ranges. The size of this dictionary was further reduced by using the low-rank basis discovered from such that more increments could have been used if required; it was found that exceeding 100 increments for each parameter was unnecessary, with changes below the noise level when trialed on the acquired data. Unless stated otherwise, all reduced dictionaries used 1Z = 1s, 2 = 84ms and 2 = 8.28μs (GM-WM average in Section 3.2).
In addition to parameter estimation, reconstructed data (̃) can be projected back to the time domain (x) from which familiar semiquantitative metrics, MTR and ihMTR, can be calculated. The degree of contrast varies through time (since the signals are time-varying), hence we define these metrics using the minimum and maximum signal values: where 0 to 3 (as indicated in Figure 1) are times within the cycle.

Imaging Experiments
All experiments were performed on a 1.5T Philips Ingenia MRI system with a 15-channel head coil for signal reception.

System Characterization and Correction for RF Instability
Our sequence was first characterized using 1D phantom scans (test-tube phantom with no phase encoding), allowing direct measurement of the full temporal signal with no need for low-rank reconstruction. Inconsistencies were found between the on-resonance amplitudes and phases of each RF pulse type, causing oscillations in signal when switching between them. The origin of these inconsistencies could not be determined but were stable for a given pulse configuration. Empirically determined complex scaling factors for each pulse were used to ensure the on-resonance component of each pulse type was the same. See Supporting Information for more details.

Sequence Implementations
For MT-MRF, a 3D tiny golden angle radial, "stack-of-stars" k-space trajectory with one spoke per RF pulse was used for all experiments (14), where TR = 5.3ms, TE = 2.7ms and FA = 29.5˚ (additional sequence details below). Encoding was ordered so as to loop over repeats of each radial 'star', before then looping over partitions. Although not necessary for the method to work, the total number of radial spokes acquired per partition was set to be commensurate with N cycle so that the same -sampling was used for all partitions, enabling independent reconstruction of each z-location. As a reference, ss-ihMT (10) was also used, in this case with Cartesian k-space sampling, acquiring separate images using 1-band, 2-band and 3-band RF pulses.

Phantom Validation
MT-MRF data were acquired on three phantoms: MnCl2-doped water (0.05mM), cross-linked bovine serum albumin (BSA, 10% w/w in water prepared as in Koenig et al. (31)) and prolipid 161 ('PL161', 15% w/w; Ashland Inc, Covington, Kentucky USA; prepared as in Swanson et al. (2) but excluding any T1 reducing agent). These phantoms are examples of "no semisolid", "no dipolar order", and "ihMT" materials respectively. FOV of 115 × 115 × 150mm and isotropic resolution 1.5mm 3 were used with a total acquisition time of 21m32s. Quantitative maps of 1 , f and 1D were computed using the reduced dictionary. In practice, it was observed that the large spread of 2 and 1 between phantoms led to inconsistent results, and so a second fit was performed using an alternative dictionary with the same fixed parameters except 2 = 130ms. Quantitative parameter estimates were compared against ss-ihMT results from the same phantoms, acquired as part of a previous study (10). That study used the same resolution but acquired data over multiple excitation flip angles with a significantly longer total scan duration of 60m48s. Quantitative maps were produced by dictionary matching using separate dictionaries relevant to this sequence.

In vivo Study
Human scanning was performed on 5 healthy subjects (3 males and 2 females, age 29±7. For comparison, ss-ihMT datasets were also acquired for every subject. These acquisitions used the same resolution, FOV, TR, TE and FA as MT-MRF scans; 6 signal averages were used to match acquisition duration (each average = 1m23s). For ss-ihMT data, motion correction using rigid body registration (imregister in MATLAB 2019a) was run before signal averaging and subsequent calculation of MTR/ihMTR as in (10). Spatial smoothing with Gaussian kernels (standard deviation 1.05mm for MTR and 1.35mm for ihMTR) was used to improve SNR. Datasets for each subject were aligned and regions-of-interest (ROIs) drawn in corticospinal tracts (CST), frontal WM, and cortical GM for numerical analysis. ss-ihMT can produce MTR/ihMTR but no other parameter estimates. Figure 2 shows ihMT contrast efficiency (η) calculated for a range of acquisition parameters.

Simulation of Obtainable Contrast
The colored areas are excluded since they violate constraints on peak B1, duty-cycle or timing. A clear region of maximum η exists between FA = 25-40˚ and ∆ = 7-9kHz. Given the optimal FA and ∆, Figure 2b shows some example sequences with different 1B and MB in a single cycle, marked by colored dots. The purple scheme has maximal contrast, though the dot marked 'scan' (yellow) was used for our experiments, with sequence parameters: TR = 5.3ms, FA = 29.5˚, τ = 2.5ms, ∆ = 8.1kHz, MB = 1B = 300 and a constant B 1,rms = 4μT (maximum achievable). This scheme was selected because it yielded more gradual contrast variations than the 'maximum contrast' case and is expected to provide favorable properties for reconstruction (24,32). A 'forbidden' scheme was also simulated (green) that exceeds maximum B1 constraints but generates even greater ihMT contrast (∆ihMT).
Simulated signal profiles are also compared in Figure 2 according to the protocols highlighted in the above heat maps, and with an ss-ihMT scan with equal time-averaged RF power (B 1,rms ). One cycle of each sequence is shown using ordering 2B-1B-3B-1B. ∆ihMT is found by comparing the end of the 2B period with the end of the 3B period ( 1 and 3 , as labelled in Figure 1). The 'scan' sequence generates almost twice the contrast achieved using ss-ihMT. The respective efficiencies of contrast generation are η = 0.080s -1 and η = 0.146s -1 . The 'maximum contrast' sequence yields only a marginally (∽6%) greater ihMT effect, though the 'forbidden' scheme contrast is significantly higher.  Figure   4b) of a dictionary corresponding to the 'scan' sequence. The first five components describe 99.2% of the 'energy' within the dictionary. Figure 4c plots the mean residuals in representing WM, GM, and CSF 'test' signals for different R. Improved representation of each tissue type occurs for higher R as expected but by R = 5, all three tissue types are accurately recreated with negligible residuals.  phantom, compared to results from a non-phase-encoded experiment (NPE); the signal is faithfully represented by the low-rank reconstruction (these data are not fitted to one another, only over-plotted). Figure 6 shows MTR and ihMTR maps of the brain from all five subjects. Table 1 and   Supporting Information Table S2 compare these metrics between each method for matched ROIs after image registration. Supporting Information Figure S4 is a test-retest comparison for the same subject, scanned eight months apart. Figure 6: Example single-slice in vivo MTR and ihMTR maps from ss-ihMT and MT-MRF and each subject (separate columns). The former shows strong GM-WM contrast, but the latter consistently appears more correlated with WM. corticospinal tracts in the axial slices from five subjects (1-5; left-to-right) in Figure 6. Bottom: MT-MRF quantitative parameter estimates obtained for the same axial regions-of-interest from the five subjects in Figure 7.   (Table 2) compares values from MT-MRF against fitting to ss-ihMT data using multiple flip angles (details in Ref. (10)). Good agreement is observed between and 1D estimates, with MT-MRF yielding tighter distributions (lower variance) even though the ss-ihMT data were acquired over almost three times longer examination. Two sets of fixed parameters were used to estimate phantom properties ( 2 = 84ms and 2 = 130ms; all others equal). Estimates of and 1D are less sensitive to this change and agree within precision when using either fixed parameter set. Though for BSA and PL161, 1 changes significantly when fixed 2 is altered. Corresponding MT-MRF parameter maps for these two scenarios are shown in Supporting Information Figure S3.  In vivo parameter maps from Subject 1 are shown in Figure 8 for single axial, coronal and sagittal slices. Figure 7 presents consistent parameter maps across all subjects, with and 1D

Quantitative Parameter Estimation
showing apparently distinct contrasts. To demonstrate this, a joint-histogram of and 1D in GM and WM pixels of Subject 1 (masks created using FSL5.0 BET and FAST (33)) is included in Figure 8. Both semiquantitative metrics ( Figure 6, Table 1, and Supporting Information Table S2) and quantitative metrics (Figure 7, Table 1) show inter-subject and test-retest (Supporting Information Figure S4) repeatability, with MT-MRF yielding a higher ihMTR across all chosen brain regions.
Supporting Information 5 presents a series of alternative fits and brief bias investigation as fixed 1Z , 2 , 2 and values are changed. Since tissues with ∽ 0 provide no MT-mediated contrast, 1 cannot be estimated for these -this is the case for CSF and MnCl2-doped water. When presenting these results, voxels in the 1 map with < 0.04 are set to zero; the threshold here was determined by examining Figure 8, which shows that all GM and WM pixels are expected to exceed this value.

DISCUSSION
This work uses a balanced SSFP sequence with multiband RF pulses that are periodically switched to allow variable off-resonant saturation (i.e. none, single frequency, and dual frequency) while keeping the on-resonance flip angle constant. MT-MRF never directly breaks the steady-state of the water magnetization but MT-mediated signal fluctuations of around 50% are observed in vivo.
These transient signal changes are reconstructed using techniques proposed for MR fingerprinting.
Indeed, as presented, the sequence is essentially a very simple MRF acquisition that would yield constant signals (and hence be unable to perform any parameter estimation) in the absence of MT.
The resulting measurements are used for quantitative estimation of some of the underlying model parameters via dictionary matching but also to reconstruct semiquantitative MT/ihMT ratio maps.

MT-MRF for Semiquantitative Mapping
Though reconstruction of artefact-free time domain images is not usually associated with classic MRF methods, the low-rank inversion (15) approach does produce high-quality 'singular' images (̃) that would normally be passed directly to dictionary matching. It is also possible to map these back to the time domain (x) and we have used this property to produce MTR/ihMTR maps using points of maximum contrast within the time-evolving signals as references (Equations 3-4).
The initial motivation for developing MT-MRF was to improve the contrast achieved from the previously developed ss-ihMT method (10). While ss-ihMT used separate image acquisitions with and without off-resonant saturation, it was shown by Varma et al. (12) that a 'low duty-cycle' scheme, alternating between high and low saturation power, would boost contrast. Sequence parameters were optimized with this in mind, trying to maximize the peak ihMT contrast (normalized to √TR). The extent of this improvement is shown in Figure 2 -ihMT contrast is almost doubled versus ss-ihMT. Figure 5 shows the reconstructed images from a phantom experiment in both the reduced basis (̃) and time domain ( ), including a comparison of reconstructed time domain signals with non-phase-encoded measurements (Figure 5d). The low-rank inversion reconstruction (using R = 5) is able to faithfully represent the time domain signal. The resulting MTR map is non-zero in PL161 and BSA phantoms, while the ihMTR map is non-zero only for PL161, as anticipated (10). Figure 6 shows MT-MRF MTR and ihMTR maps from a study on five healthy subjects, compared to those from ss-ihMT. Contrast is consistent across subjects (Tables 1 and S2), whilst a test-retest validation on one subject (Supporting Information Figure S4; Table S3) also showed very consistent contrast. MT-MRF generally finds greater ihMTR than for ss-ihMT (62% greater in frontal WM, for example) but lower MTR; since the scheme was optimized for ihMT, the signal never fully returns to equilibrium. Standard deviations in Table 1 are smaller for MT-MRF. This may be due to a combination of better CNR, superior motion robustness since a radial readout was used, and use of spatial regularization (Equation 2).

Quantitative Parameter Estimation
The generation of MTR and ihMTR maps is a limited use of the time-series data available; quantitative parameter estimation is also possible. CRLB calculation (Figure 3) suggested that not all parameters can be estimated to good precision, thus it was decided to only reconstruct 1 , and 1D .
Phantom experiments using MT-MRF were compared to previous measurements using ss-ihMT with variable flip angles for estimation of these three parameters ( Table 2). Estimates of f and 1D agree for both methods, and the MT-MRF sequence provides better precision, particularly for 1D . The estimated 1D ~ 25ms for PL161 is in line with measurements from Swanson et al (2). Fits using different 2 were performed because BSA and PL161 have very different properties (more so than brain tissue) and thus fixing parameters was found to cause different 1 estimation bias for MT-MRF and ss-ihMT. However, MT parameters estimates (i.e. f and 1D ) were consistent regardless of 2 .
In vivo results in Figure 7 show consistent parameter values across subjects. and 1D are significantly different in GM compared to WM and 1D seems to highlight strongly myelinated structures (e.g. CST). Generally, parameter estimates between different subjects are in agreement (Table 1), with the highest 1D values found in CST (~5.5ms) and lower values (~2.2ms) found in cortical GM; orientation with respect to B 0 may influence these values (34). Estimates for exhibit higher precision -an observation that is supported by the CRLB analysis in Figure 3. Maps of and 1D do show distinct contrasts; Figure 8 presents a whole-brain joint-histogram from one subject and demonstrates that these parameters are correlated but not simply linearly dependent.
Compared to the literature, CST 1D estimates (∽6.0ms in the coronal slice in Figure 8) are generally in agreement with those reported elsewhere (35) and WM estimates for concur with those in Varma et al. (12). Others (35,36) have reported GM 1D to be similar to WM (∽5.6ms), however Swanson et al. also found far shorter 1D in GM compared to WM from bovine spinal cord. (2)

Parameter Estimation Biases
Parameter fixing is common in quantitative MT imaging because model parameters can influence the signals in similar ways, so estimation is not well-posed and may create bias. In our work, and 1D estimates are comparable to those found by others. The phantom experiment provided similar estimates from two different experiments, though separate fixed values of 2 were considered. The main impact of changing 2 is to alter the estimate of 1 , which could indicate that this is acting as a 'nuisance' parameter to absorb uncertainties in the relaxation times.
Supporting Information Figure S5c explores this further for in vivo data by examining the effects of changing each fixed model parameter by ±10%. It is observed that 1Z has negligible impact, while changing causes small changes in and 1D . In line with phantom experiments, changing 2 causes different 1 estimates but has negligible impact on and 1D . The most significant correlation exists between 1D and 2 , which has also been reported by Varma et al. (35).
Recent work from Wang et al. (39) has questioned the commonly made assumptions that 1Z can either be fixed at some simple value (i.e. 1s) or set equal to 1 (40,41). Instead, Wang et al. suggest that 1Z is a major determinant of observed T1 and is strongly frequency dependent; a 1Z of approximately 120ms in WM should be expected at 1.5T. Figure 8 used the assumption 1Z = 1s but we repeated in vivo dictionary fitting for Subject 1 using a fixed value of 200ms; the result in Supporting Information Figure S5a shows substantially longer 1 but largely unchanged and 1D . CRLB analysis (Figure 3) showed that 2 could be estimated instead of 1 -this was also performed (fixing 1 = 1.36s from Section 3.2) and the result (Supporting Information Figure S5b) shows that a plausible T2 map can be obtained but again, there is little impact on estimated and 1D .

Determining Subspace Rank
The embodiment of MT-MRF used in this work was chosen to emulate a 'standard' MT measurement that cycles over different saturation settings with the same flip angle; it has the interesting property that all signal variation is MT-mediated. One important side-effect of this is that parameter estimation is not possible for substances without a semisolid compartment. However, the fairly simple time-structure of this specific sequence leads to smooth temporal variations such that the generated dictionaries have an effective rank of five. Time domain signals from representative tissues could be reconstructed with <0.75% error using the reduced basis. The presented images used rather long acquisitions (~25 minutes) which were necessary because ihMT contrast is rather weak.
In one sense, the images are highly undersampled since approximately nine complete volumes of k-space data were acquired in order to reconstruct a time-variable cycle of 1200. But given the effective rank of five, the reconstruction problem is actually over-determined, though the effective degree of undersampling does not map uniformly across the singular images.
By its nature, the sequence does not probe all degrees of freedom available. This may explain the low rank and also the inability to measure more parameters to high precision. Nevertheless, we were able to measure some important parameters such as f and 1D . Estimation accuracy is hard to assess because no 'gold standard' measures exist (particularly for ihMT) but our findings are similar to those reported elsewhere. Though 1 was also estimated, analysis showed that this parameter may be more biased due to the fixing of other relaxation times in the current implementation. Beyond this proof-of-concept demonstration, future work will consider combining the multiband pulses used here with more normal MRF sequences that also modulate flip angle and include other contrast preparation pulses. Modulating RF frequency offsets may allow us to disentangle the interdependence between 1D and 2 . These sequences could be optimized for precision (rather than maximum ihMT contrast) using established CRLB based methods (28) to boost parameter estimation capability, and a more extensive analysis could be performed to identify and reduce remaining estimation biases.

CONCLUSIONS
This work presents a cyclic-steady-state sequence using multiband pulses with modulated off-resonance power. As presented, the method is a proof-of-concept that uses the machinery of modern MR fingerprinting on a sequence with constant flip angle, that should not be able to perform any parameter estimation if the normal Bloch Equations model is used. We show that MT-mediated signal changes are actually large and can be used to permit semiquantitative and limited quantitative parameter estimation. Future embodiments of this method will involve more general integration with MRF, using the multiband pulse off-resonance bands as an additional parameter to vary alongside others already used for fingerprinting.

Data Availability Statement:
The code used to generate the simulation results in this study can be found at: https://github.com/mriphysics/MT-MRF (hash 068cbdac was used for the presented results).

Supporting Information 1: Steady-State Signal Calculation
The time-evolution of magnetization can first be written as Eq. S1 and then reformulated as Eq. S2,

Supporting Information 2: Multiband Pulse Corrections
MT-MRF switches between single-band and multiband pulses with the aim that on-resonance components of these pulses are invariant. When conducting non-phase-encoded experiments, we noticed that significant and unexpected signal fluctuations occurred at the point where pulse type changed and found this was caused by the on-resonance lobes having a slightly different amplitude and phase. A pick-up coil measurement utility already implemented on the scanner was used to characterize these discrepancies for each scan. The origin could not be found but was suspected to be from the RF hardware. Errors were not easily predictable from scan parameters but were found to be consistent as long as scan parameters were unchanged.
Example measurements are plotted in Figure S1 that displays the measured RF pulses in the frequency domain. When zooming in on the on-resonance band, it is clear that there are some discrepancies between the pulses. Complex scaling factors were calculated to make the behavior of each multiband pulse (i.e. 2B and 3B) match that of the single-band (1B) pulse at Δ = 0. The MT-MRF pulse sequence was implemented to allow user-defined values to be entered for these scaling factors. Table S1 shows examples of the required scaling factors. Figure S1: Plots used to derive amplitude and phase correction factors between single-band (SB), 2-band (2B) and 3-band (3B) pulses. Note how the on-resonance lobes and phases of each pulse are slightly offset from one another. These offsets are overcome using the example factors shown in Table S1. Non-phase-encoded phantom scans were then repeated to ensure that fluctuations were removed, and this is confirmed in Figure S2 (below). Figure S2: Example 1D signal profiles from each phantom before (unscaled, US) and after (scaled, S) amplitude and phase correction. Large fluctuations at transitions between pulse types are significantly reduced and removed for BSA and PL161. The pulse cycle used to obtain these plots is: {500 2B + , 500 SB, 500 2B -, 500 SB, 500 3B, 500 SB etc.} so does not match the actual acquisition scheme used in the manuscript; it was used for debugging purposes only. 2B + refers to a 2-band pulse with a positive off-resonance lobe and 2Bhas an equal but opposite negative off-resonance lobe.
Though the latter is not used in our phantom and in vivo experiments, it can be used to generate ihMT contrast when combined with a matched 2B + pulse to approximate dual frequency off-resonance.

Supporting Information 3: Supplementary Phantom Results
To supplement Table 2 in the manuscript, below are single-slice maps resulting from the two separate dictionary fits to MT-MRF phantom data. Note the similarity between MT parameter estimates when either 2 = 84ms or 2 = 130ms. Only 1 values change slightly due to an expected 2 dependence. Figure S3: Top: MT-MRF parameter maps obtained using a GM-WM average parameter set with 2 = 84ms (Section 3.2). Bottom: Equivalent parameter maps but assuming a more 'PL161-like' fixed value for 2 (130ms) only.

Supporting Information 4: Supplementary in vivo Results
To supplement the inter-subject ihMTR comparison shown in Table 1 of the manuscript, below is an equivalent comparison of corresponding MTR values between MT-MRF and ss-ihMT. Table S2: MTR values obtained from different regions-of-interest for the axial slices of the five subjects in Figure 6.

PRE-PRINT: SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE
For Subject 2 (male, aged 25), two identical MT-MRF scans were completed approximately eight months apart. Figure S4 presents a central axial slice from each dataset following registration. Table   S3 compares measurements in corresponding matched regions-of-interest. Figure S4: Equivalent central axial slices from the same healthy subject but two separate acquisitions. ihMTR and MTR contrasts from MT-MRF are repeatable according to the region-of-interest measurements shown in Table S3. Table S3: ihMTR and MTR values from three different brain regions in the axial slices shown in Figure S4.

Fixing 1Z to a Lower Value
Recently, Wang et al. suggested that 1 is much lower than is usually assumed in MT literature. (2) Therefore, we repeat dictionary fits for Subject 1 using a lower, brain-average 1Z = 200ms. Figure S5a: Dictionary fitting results for Subject 1 assuming lower 1Z . Compared to Figure 8, 1 estimates increase since free and semisolid T1 are coupled to one another, whilst estimates for and 1D are mostly unchanged.

Fitting for 2
From Figure 3, Combination 4 seems to also give reasonable estimation precision and so dictionary fits are repeated for Subject 1 but for the case where 2 is estimated and 1 is fixed at 1.36s. (3) Figure S5b: Dictionary fitting results for Subject 1 assuming fixed 1 but estimating 2 . 2 maps show similar contrast to those for 1 but once more, estimates for and 1D are very similar to those reported in Figure 8.

PRE-PRINT: SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE
Dictionary fitting was also performed by incrementing our previously assumed fixed parameter values by ±10%. Difference maps were then computed between these two to reveal possible interdependence of model parameters and investigate the sources of bias in our in vivo dictionary fits. Figure S5c: Difference maps for each estimated parameter given a ±10% change in the fixed parameters stated on the left-hand side. Each estimated parameter seems invariant to small changes in 1Z ; 1 is most coupled to 2 ; 1D shows considerable dependence on 2 ; and is slightly influenced by . Colorbars for each estimated parameter change (∆X, where X is the estimated parameter) are chosen to cover 60% of the ranges used for these quantities in other figures.