Unconstrained quantitative magnetization transfer imaging: disentangling T 1 of the free and semi-solid spin pools

Since the inception of magnetization transfer (MT) imaging, it has been widely assumed that Henkelman’s two spin pools have similar longitudinal relaxation times, which motivated many researchers to constrain them to each other. However, several recent publications reported a T s 1 of the semi-solid spin pool that is much shorter than T f 1 of the free pool . While these studies tailored experiments for robust proofs-of-concept, we here aim to quantify the disentangled relaxation processes on a voxel-by-voxel basis in a clinical imaging setting, i.e., with an effective resolution of 1.24mm isotropic and full brain coverage in 12min. To this end, we optimized a hybrid-state pulse sequence for mapping the parameters of an unconstrained MT model. We scanned four people with relapsing-remitting multiple sclerosis (MS) and four healthy controls with this pulse sequence and estimated T f 1 « 1 . 84s and T s 1 « 0 . 34s in healthy white matter. Our results confirm the reports that T s 1 ! T f 1 and we argue that this finding identifies MT as an inherent driver of longitudinal relaxation in brain tissue. Moreover, we estimated a fractional size of the semi-solid spin pool of m s 0 « 0 . 212, which is larger than previously assumed. An analysis of T f 1 in normal-appearing white matter revealed statistically significant differences between individuals with MS and controls.


Introduction
Longitudinal relaxation is a vital contrast mechanism in magnetic resonance imaging (MRI).For example, the MP-RAGE (Mugler and Brookeman, 1990) pulse sequence generates excellent gray matter (GM)-white matter (WM) contrast and-compared to mostly T 2 -weighted pulse sequences like FLAIR (Hajnal et al., 1992)-may be more specific to the underlying tissue changes in multiple sclerosis (MS) lesions (Barkhof, 1999;Bagnato et al., 2003).Koenig et al. (1990) discovered that macromolecules and lipids, in particular myelin, are the source of fast longitudinal relaxation in WM.Though their experiments were not designed to identify the mechanism through which macromolecules facilitate relaxation, they hypothesized that magnetization transfer (MT) (Wolff and Balaban, 1989) is a driving force of relaxation.
Magnetization transfer is commonly described by Henkelman's two-pool model (Henkelman et al., 1993), where one spin pool, the free pool, consists of all protons bound in water and is denoted by the superscript f , and the other pool, the semi-solid pool, consists of protons bound in macromolecules (e.g., proteins and lipids) and is denoted by the superscript s.In standard clinical pulse sequences, one does not directly observe the latter spins because their transversal magnetization relaxes below the noise level before it can be observed (T s 2 « 10µs).However, the exchange of longitudinal magnetization between the two pools alters the free pool's longitudinal magnetization, resulting in bi-exponential relaxation.
The indirect nature of the semi-solid spin pool's impact on the MRI signal entails an entanglement of differ-ent parameters.In order to mitigate the consequent noise amplification, most quantitative MT (qMT) approaches constrain T s 1 " T f 1 resulting in T s 1 « 1.1s (Yarnykh, 2002;Dortch et al., 2011) or, similarly, simply assume that T s 1 " 1s (Henkelman et al., 1993;Morrison and Henkelman, 1995).However, recent studies have suggested that T s 1 « 0.3s and T f 1 « 2s in white matter at 3T (Helms and Hagberg, 2009;Gelderen et al., 2016;Manning et al., 2021;Samsonov and Field, 2021).These studies overcame the challenges of an unconstrained model by either using brain-wide estimates of T s 1 and/or T f 1 (Gelderen et al., 2016;Samsonov and Field, 2021) or fitting the MT model to NMR samples (Manning et al., 2021) or a single large ROI averaged over multiple participants (Helms and Hagberg, 2009).
Our work aims to confirm these findings and to offer evidence in support of Koenig's hypothesis that MT is a key driver of longitudinal relaxation in brain tissue.Moreover, we aim to provide, for the first time, voxel-wise fits with the unconstrained two-pool MT model.Key to this advance is a hybrid state (Assländer et al., 2019b) of the free spin pool that can provide increased efficiency in the en-coding and the disentanglement of the MT and relaxation processes (Assländer, 2021).Further, we describe the semi-solid spin pool with the generalized Bloch model for slight improvements in model accuracy (Assländer et al., 2022).
We first validated the approach with phantom scans.Then, we measured reference parameters in vivo using 36min scans in participants with multiple sclerosis and healthy controls.Last, we tested rapid imaging protocols and found that our proposed approach enables unconstrained qMT imaging with an effective resolution of 1.24mm, 1.6mm, and 2.0mm isotropic in 12, 6, and 4 minutes, respectively.

Magnetization Transfer Model
We use the MT model described in Assländer et al. (2022Assländer et al. ( , 2024)), which builds on Henkelman's two-pool spin model (Henkelman et al., 1993) and captures the two pools with a Bloch-McConnell equation (McConnell, 1958): The free pool, sketched in red in Fig. 1, captures all protons bound in liquids where fast molecular motion causes an exponential relaxation of the transversal magnetization with a characteristic T f 2 Á 50ms (Bloembergen et al., 1948).The free pool's magnetization is described by the Cartesian coordinates x f , y f , z f , the off-resonance frequency is described by ω z , and the Rabi frequency of the RF pulses by ω y .For readability, we here use relaxation rates (R f,s 1,2 " 1{T f,s 1,2 ).The magnetization components x s , z s of the semi-solid spin pool, sketched in purple in Fig. 1, capture all protons bound in large molecules such as lipids.The motion of such molecules is restricted, resulting in a much faster and non-exponential relaxation with a characteristic time constant of T s 2 « 10µs, which prevents a direct detection of this pool with standard clinical MRI.Within the brain parenchyma, we assume the decay characteristics associated with a super-Lorentzian lineshape (Morrison and Henkelman, 1995).The nonexponential characteristics of this lineshape prohibit a description with the original Bloch equations, but such dynamics can be described with the generalized Bloch model (Assländer et al., 2022).For computational efficiency, we can approximate the non-exponential decay by an effective exponential decay with a linearized relaxation rate R s,l 2 pR s 2 , α, T RF q.While exponential and non-exponential decays necessarily deviate, we can identify an R s,l 2 that results in the same magnetization at the end of an RF pulse.To this end, R s,l 2 depends on the flip angle α and the duration T RF of respective RF-pulse in addition to the biophysical parameter R s 2 .We neglect the y s component assuming, without loss of generality, ω x " 0 and given that R s,l 2 " ω z .The exchange rate R x captures exchange processes between the pools.A sixth dimension is added to allow for a compact notation of the longitudinal relaxation to a non-zero thermal equilibrium.
Throughout the literature, multiple normalizations of m s 0 have been used.Here, we use m f 0 `ms 0 " 1 so that m s 0 describes the fraction of the overall spin pool, a definition which has also been used, e.g., by Yarnykh (2012).Other papers, such as Henkelman et al. (1993); Gochberg and Gore (2003), measure ms 0 " m s 0 {m f 0 or, equivalently, normalize to mf 0 " 1.The conversion between the two definitions is simply m s 0 " ms 0 {p1 `m s 0 q and ms 0 " m s 0 {p1 ´ms 0 q.
2 free / liquid Figure 1: Sketch of the two-pool magnetization transfer model (Henkelman et al., 1993).This model jointly describes all magnetization arising from protons bound in liquids by the spin pool m f 0 , and all magnetization arising from protons bound in macromolecules by the pool m s 0 whose transversal relaxation time is several orders of magnitude shorter.We normalize the thermal equilibrium magnetization to m f 0 `ms 0 " 1 and describe the magnetization transfer between the pools by the rate Rx " 1{Tx.The model is governed by Eq. (1).

Comparison to constrained MT models
In the absence of RF pulses (ω y " 0), we can isolate the longitudinal components of Eq. ( 1): (2) An eigendecomposition of the Hamiltonian in Eq. ( 2) has three distinct eigenvalues (Henkelman et al., 1993;Gochberg and Gore, 2003;Yarnykh, 2012).One is zero and corresponds to thermal equilibrium.The smaller remaining eigenvalue (in the absolute value) can be considered an apparent relaxation rate of the free pool R f,a 1 that is approximated by the following Taylor expansion at R s The MT contributions to R f,a 1 therefore depend foremost on the macromolecular pool size m s 0 and the two relaxation rates.Higher order terms additionally depend on the exchange rate R x .Eq. ( 3 1 , making MT an important driver of longitudinal relaxation.For example, let us assume R f 1 " 0.5{s, R s 1 " 3/s, and m s 0 " 0.2.In this case, the linear correction term is 0.5/s and, thus, R f,a 1 « 1.0{s ff R f 1 .The largest eigenvalue (in absolute value) is given by which is dominated by the exchange rate R x for many tissues.Hence, it can be interpreted as a cross-relaxation term and we henceforth refer to it as the apparent exchange rate.
From the eigenvectors, we can derive a Taylor expansion of the apparent semi-solid pool size (Appendix A): (5) Eq. ( 5) reveals that m s 0 is underestimated when assuming R s 1 " R f 1 .To give a sense of the magnitude of this bias, we can insert the above example values and further assume R x " 15{s, which results in m s,a 0 « 0.15 instead of the underlying m s 0 " 0.2.

Pulse sequence design
As mentioned above, we utilize the hybrid state (Assländer et al., 2019b) and its flexibility to encode and disentangle the different relaxation mechanisms.Similar to balanced SSFP sequences (Carr, 1958), we balance all gradient moments in each T R .On the other hand, we vary the flip angle and the duration of the RF pulses.During slow flip angle variations, the direction of the magnetization establishes a steady state and adiabatically transitions between the steady states associated with different flip angles.As we showed in Assländer et al. (2019b), moderate change rates of the flip angle simultaneously yield a transient state of the magnetization's magnitude, and we call this combination the hybrid state.It combines the tractable off-resonance characteristics of the bSSFP sequence, particularly the refocusing of intra-voxel dephasing (Carr, 1958;Scheffler and Hennig, 2003), with the encoding potential of the transient state.
Our pulse sequence consists of a rectangular π inversion pulse, flanked by crusher gradients, followed by a train of rectangular RF pulses with varying flip angles and pulse durations.The RF phase is incremented by π between consecutive RF pulses.The pulses are separated by a T R " 3.5ms, which is approximately the minimal T R with which we can perform gradient encoding with |k max | " π{1mm and avoid stimulating the peripheral nerves.After 1142 RF-pulses, i.e., after a cycle time of 4s, the remaining magnetization is inverted by the subsequent π pulse, then the same pulse train is repeated.
The relaxation and MT processes are encoded with two established mechanisms.First, the T 2 -selective inversion pulse inverts the free pool while keeping the semi-solid pool largely unaffected.As described by Gochberg and Gore (2003), this induces a bi-exponential inversion recovery curve of the free pool composed of its intrinsic longitudinal relaxation and cross-relaxation to the semi-solid spin pool.Second, we can use the flip angle and the pulse duration to control the different relaxation paths.In good approximation, the RF-pulse duration only affects the saturation of the semi-solid spin pool's longitudinal magnetization (Gloor et al., 2008).In contrast, changes in the flip angle affect the relaxation processes of the free pool (Assländer et al., 2019a,b), the magnetization transfer between the two pools, and the saturation of the semi-solid spin pool (Gloor et al., 2008).More details on this interplay can be found in Assländer et al. (2024).

Numerical optimization of the pulse train
We numerically optimized the flip angles and pulse durations of RF-pulse trains based on these two encoding mechanisms.The optimization objective was the Cramér-Rao bound (CRB) (Rao, 1945;Cramér, 1946), which predicts the noise variance of a fully efficient unbiased estimator.We note that least squares fitting and the neural network-based fitting used in this article (cf.Section 3.6) are, strictly speaking, neither fully efficient nor unbiased (Newey and McFadden, 1994;Wu, 1981).Nonetheless, the CRB can be used as a proxy for the "SNR-efficiency" or "conditioning" (Zhao et al., 2019;Haldar and Kim, 2019) and we adopt this heuristic here.
We calculated the CRB as described in Assländer et al. (2024) and optimized for the CRBs of the relaxation rates and the other model parameters.We optimized a separate pulse train for each of the biophysical parameters 1 , and T s 2 , while additionally accounting for ω z , B 1 " ω y {ω nominal y , and the scaling factor M 0 as unknowns, where M 0 jointly describes the overall spin density and receive-coil sensitivity profiles.Additionally, we optimized a pulse train for the sum of the CRBs of all biophysical parameters, normalized with respective squared parameter values to resemble the inverse squared SNR.We performed all simulations and CRB calculations with m s 0 " 0.25, R f 1 " 0.5{s, R f 2 " 15.4{s, R x " 20{s, R s 1 " 2/s, T s 2 " 10µs, ω z " 0, and B 1 " 1.The resulting spin trajectories (Fig. S1) and the corresponding CRB values (Tab.S1) are discussed in the supporting information Section S1.Supporting Section S2 and Figs.S2, S3, and S6 connect the CRB to experimental noise measurements.

Phantom scan
We built a custom phantom composed of cylindrical 50mL tubes filled with different concentrations of thermally cross-linked bovine serum albumin (BSA).We mixed the BSA powder (5%, 10%, . . ., 35% of the total weight) with distilled water and stirred it at 30 ˝C until the BSA was fully dissolved.We filled 7 tubes with the resulting solutions and thermally cross-linked them in a water bath at approximately 90 ˝C for 10 minutes.
We scanned this phantom on a 1.5T Sola and 2.9T Prisma scanner (Siemens, Erlangen, Germany).We performed a 6min scan with each of 6 individual optimizations, resulting in a 36min overall scan time.For each 6min scan, the RF pattern is repeated 90 times, during which we acquire 3D radial k-space spokes with nominal 1.0mm isotropic resolution (defined by |k max | " π{1mm).The sampled k-space covers the insphere of the typically acquired 1.0mm k-space cube.By comparing the covered k-space volume, we estimate an effective resolution of 1.24mm, which we report throughout this paper (Pipe et al., 2011).We note that the stated effective resolution does not account for blurring introduced by undersampling in combination with a regulized reconstruction.We changed the direction of the k-space spokes with a 2D golden means pattern (Winkelmann et al., 2007;Chan et al., 2009) that is reshuffled to improve the k-space coverage for each time point and to minimize eddy current artifacts (Flassbeck and Assländer, 2023).

In vivo scans
Each participant's informed consent was obtained before the scan following a protocol that was approved by the NYU School of Medicine Institutional Review Board.To establish high-quality reference data, we performed in vivo scans of 4 individuals with clinically established relapsingremitting MS (extended disability status scale (EDSS) 1.0-2.5, unknown for one participant; no recent history of relapses; age 37.5 ˘8.7; 3 female) and 4 healthy controls (age 28.8 ˘5.6; 3 female) with a 2.9T Prisma scanner and the 36min protocol described in Section 3.3.In addition to the hybrid-state scans, we acquired 3D MP-RAGE and FLAIR scans, each with a 1.0mm isotropic resolution.
To test more clinically feasible scan times, we scanned an additional participant with MS with three "rapid" protocols with different effective resolutions:

Image reconstruction
We performed retrospective motion correction similar to Kurzawski et al. (2020).However, our approach deviates from Kurzawski et al. in one key aspect: Instead of using an SVD to maximize the first coefficient's signal intensity, we utilize a generalized eigendecomposition (Kim et al., 2021) to maximize the contrast between brain parenchyma and CSF (Flassbeck et al., 2024).We reconstructed images directly in the space spanned by three basis functions associated with the generalized eigendecomposition (Tamir et al., 2017) and used a total variation penalty along time to reduce undersampling artifacts (Feng et al., 2014).The reconstructions were performed with a spatial resolution of 4mm isotropic and a temporal resolution of 4s.The images corresponding to the first coefficient were co-registered using SPM12 and the extracted transformation matrices were subsequently applied to the k-space data (translations) and trajectory (rotations) to correct the full-resolution reconstruction.
For the 36min in vivo scans, we used separate subspaces for each 6min sub-scan, implemented as a blockdiagonal matrix to permit joint regularization.For the phantom scan and the rapid protocols, we reconstructed all data of the 6 sub-scans into a joint 15-dimensional subspace with otherwise identical settings.

Model fitting
For computational efficiency and robustness, we used neural networks to fit the MT model, voxel by voxel, to the reconstructed coefficient images (Cohen et al., 2018;Nataraj et al., 2018;Duchemin et al., 2020;Zhang et al., 2022;Mao et al., 2023a).This approach includes a datadriven B 0 and B 1 correction as detailed in Assländer et al. (2024).For each voxel, the complex-valued coefficients are normalized by the first coefficient and split into real and imaginary parts, which are the inputs to the network.The network retains a similar overall architecture to the design described in Fig. 2 of Zhang et al. (2022): 11 fully connected layers with skip connections and batch normalization, and a maximum layer width of 1024.The network estimates all 6 biophysical parameters of the unconstrained MT model.For both reconstruction protocols, we trained networks using the Rectified ADAM optimizer (Liu et al., 2019) to convergence with individually-tuned learning rates.For more details, we refer to Zhang et al. (2022) and Mao et al. (2023a).

Region of interest analysis
For the 36min reference scans, we registered the skullstripped (Hoopes et al., 2022) qMT maps and the FLAIR images to the MP-RAGE with the FreeSurfer package ("mri robust register") (Reuter et al., 2010).We also used FreeSurfer ("recon-all") to segment the brain based on the MP-RAGE and the FLAIR (Fischl et al., 2002(Fischl et al., , 2004)).We extracted region of interest (ROI) masks for the entire normal-appearing white matter (NAWM), several WM subregions, the cortical GM, and subcortical GM structures.To ensure that MS lesions were excluded from the ROIs, we calculated lesion masks with an in-house developed deep learning model based on the nnUNet framework (Isensee et al., 2021) 2. The fitted function has the form a 0 `aBSA cBSA `aB 0 ¨B0 `a2 ¨B0 ¨cBSA .The white background identifies coefficients that differ from zero at a 95% confidence interval.The gray background identifies coefficients for which such determination cannot be made.The rightmost column denotes the coefficient of determination.
to reduce partial volume effects with other tissues and to ensure that all ROI voxels are at least one voxel away from any lesion.

Phantom scan
The phantom validation aims to identify the parameters' dependency on the sample's protein (BSA) concentration and the magnetic field strength and to compare these findings to previous work and our general understanding of relaxation.To this end, we performed a generalized linear model (GLM) fit of the data with the BSA concentration (c BSA ) and the magnetic field strength (B 0 ) as independent variables (Fig. 2 and Tab. 1).
The estimates of the semi-solid spin pool size are consistent with the linear model m s 0 " a BSA ¨cBSA , i.e., we can reject neither the null hypothesis that m s 0 is independent of B 0 nor that it vanishes at c BSA " 0.
For R f 1 , our data supports the linear model R f 1 " a 0 `aBSA ¨cBSA `a2 ¨cBSA ¨B0 .I.e., the data suggests a linear dependence of R f 1 on the BSA concentration where the slope also depends on the field strength.By contrast, we did not observe a dependency of R f 1 of pure water (c BSA " 0) on B 0 .This finding is consistent with the very small change predicted by the Bloembergen-Purcell-Pound (BPP) theory (approximately 0.004% for pure water with a correlation time τ c " 5ps; (Bloembergen et al., 1948)).Further, the intercept a 0 " 0.28{s with a 95% confidence interval r0.15, 0.41s{s agrees with the 0.255{s predicted by the BPP theory for 1.5T and 2.89T (the BPP theory predicts an identical rate at both field strengths within the indicated precision).
For R f 2 , our data is consistent with a linear dependence on the BSA concentration and no dependence on B 0 (R f 2 " a BSA ¨cBSA ).The latter is also consistent with the very small change predicted by the BPP theory (approximately 0.001%).While the negative intercept a 0 " ´2.9{s is not physical, the 95% confidence interval r´7.6, 1.8s{s includes the R f 2 « 0.255{s predicted by the BPP theory.We detected no dependence of R s 1 on c BSA , but a statistically significant dependence on B 0 (R s 1 " a 0 `aB0 ¨B0 ), which is consistent with the reports of Wang et al. (2020).For R x , we observe neither a c BSA nor a B 0 dependency, and for T s 2 , the B 0 dependency is just above the 95% significance threshold.
Beyond the linear model, we observe increased variability of R x , R s 1 , and T s 2 estimates with decreasing c BSA , which likely stems from the smaller spin pool size.For this reason, we excluded the most extreme case (c BSA " 0.05) from the GLM fits as indicated by the brackets in Fig. 2.

In vivo reference scans
Fig. 3 demonstrates the feasibility of unconstrained qMT imaging with a hybrid-state pulse sequence, i.e., encoding all 6 biophysical parameters on a voxel-by-voxel basis.By comparing the qMT maps to the routine clinical contrasts, we observe overall good image quality in m s 0 , R f 1 , and R f 2 .However, the cerebellum reveals a slightly reduced effective resolution compared to the nominally equivalent resolution of the MP-RAGE (Fig. 3d vs. f,h,. . .).The R x and R s 1 maps exhibit reduced image quality, consistent with their higher CRB values (Tab.S1).Also consistent with its large CRB, the T s 2 map has the highest noise and artifact levels, which might also be, in part, due to a residual B 1 artifact caused by incomplete spoiling of the inversion pulse.We also find subtle residual B 1 artifacts in R f 2 (Fig. 3i,j, and Fig. S5c) and residual B 0 artifacts in a few voxels at the center of the bSSFP banding artifact (Fig. 3f,h,. . .at the base of the frontal cortex).Overall, however, we observe good performance of the data-driven B 0 and B 1 correction.
Among all qMT parameters, we observe the largest quantitative GM-WM contrast in m s 0 , followed by R f 1 .In R f 2 , however, we observe only a subtle contrast between cortical GM and WM.An ROI analysis confirms this finding: we estimated T f 2 " p83 ˘15qms and p76.9 ˘8.3qms for cortical GM and WM, respectively, which is a smaller difference compared to the difference between previously reported values (p99 ˘7q vs. p69 ˘3qms; (Stanisz et al., 2005)).Consistent with previous reports, we observe the shortest T f 2 " p59.3 ˘5.6qms in the pallidum (Fig. 3i).The exchange rate R x , R s 1 , and T s 2 exhibit little GM-WM contrast and we note that the most prominent contrast in R x and R s 1 occurs in voxels subject to partial volume effects and in CSF.In voxels with partial volume, the model might be inaccurate and the small m s 0 makes estimates of semi-solid spin-pool characteristics unreliable, in particular with an unconstrained model, as previously demonstrated by Dortch et al. (2018).Estimates of the unconstrained MT model's parameters are reported in Tab. 2 for selected WM and GM structures.compares the apparent qMT parameters of those fitted with the unconstrained model, with more brain ROIs analyzed in Supporting Tab.S2.Below, we discuss the salient differences in white and cortical gray matter.
The exchange rate estimated with the unconstrained MT model, R x " p13.6 ˘1.1q/s, is slightly lower compared to the corresponding literature ((18.1˘3.6)/s(Helms and Hagberg, 2009)) as is the the apparent exchange rate (Eq.( 4)): R a x " p16.1 ˘1.2q/s compared to (23 ˘4)/s (Stanisz et al., 2005).Notwithstanding, our analysis matches previous findings that the constraint T s 1 " T f 1 biases R x to larger values (Helms and Hagberg, 2009).Supporting Figs.S4, S5, and S9 compare the unconstrained qMT estimates to constrained fits of the same hybrid-state data.Most constrained estimates match neither the unconstrained estimates nor literature values.3)-( 5)) to unconstrained parameter estimates.The box plots pool all normal-appearing white matter voxels of each participant.The markers ˚and ˚˚indicate statistically significant differences at the p ă 0.05 and p ă 0.01 levels in a comparison of each subject's median qMT parameter estimates between participants with MS and controls.
Gray matter.An ROI analysis of the cortical gray matter, averaged over all healthy volunteers, reveals trends similar to the WM analysis: T f 1 " p2.46 ˘0.56qs and T s 1 " p0.42 0.40qs differ substantially from one another.The apparent T f,a 1 " p1.62 ˘0.23qs is in line with the mono-exponential estimate (1.82 ˘0.11)s measured by Stanisz et al. (2005).
As expected, the estimated m s 0 " 0.098 ˘0.026 is both smaller than that measured in WM and similar to the literature value m s 0 " 0.086 (derived from ms 0 " 0.094) estimated with the unconstrained MT model (Helms et al., 2004).The estimated m s,a 0 " 0.071 ˘0.051 is in line with literature values based on a constrained MT model (0.050 ˘0.005 (Stanisz et al., 2005)), though noise amplification resulting from Eq. ( 5) limits the value of this comparison.
The estimated R x " p14.0 ˘3.1q/s as well as R a x " p16.4 ˘3.4q/s of GM are, similarly to WM, lower than literature values that are based on a constrained MT model ((40 ˘1)/s (Stanisz et al., 2005)).and p ă 0.01 levels.We note that the panels for m s 0 , R f 1 , and Rx are a repetition of the ones in Fig. 5.
In NAWM, the median T f 2 of each MS subject averaged over all participants with MS was 2.1ms larger than in controls (p ă 0.05).When analyzing all unconstrained qMT parameters for the ROIs listed in Tab. 2, we found statistically significant changes of With the constrained MT model, we only found significant differences in T f,a 1 of the putamen (p ă 0.05).Fig. S10  depicts T f 1 for the ROIs above.

MS lesions
In MS lesions, we observe a substantial reduction of m s 0 (Figs. 8 and S5) and m s,a 0 (Fig. S7) relative to the NAWM, consistent with the expected demyelination.
When jointly analyzing T f,a 1 and m s,a 0 across all MS lesions using principal component analysis (Fig. 7), we find  that the first component explains 93% of the variability.By comparison, only 79% of variability is explained by the first component for the unconstrained model, suggesting an increase in independent information across the qMT parameters.This might be beneficial in understanding the various biophysical processes contributing to disease, which we elaborate on in the Discussion.

Rapid qMT imaging
All data described thus far were acquired with 1.24mm isotropic resolution and 36min scan time.To gauge the potential of our qMT approach for more clinically feasible scan times, we scanned an individual with MS with different resolutions and scan times.With 1.24mm isotropic resolution and 12min scan time, we observe overall good image quality despite slightly increased blurring and noise compared to the 36min scan (cf. the cerebellum in Fig. 3f to the one in Fig. 8a).With 1.6mm isotropic nominal resolution and 6min scan time, we observe similar image quality besides the reduced resolution, and the same is true for 2.0mm isotropic in 4min.

Unlocking unconstrained qMT imaging
Many pulse sequences are sensitive to R f 1 and R s 1 .Using SIR (Dortch et al., 2018) as an example, this can be illustrated with normalized CRB values (cf.supporting Tab.S1) under the assumption that only the respective parameter is unknown: norm.CRBpR f 1 q « 25s and norm.CRBpR s 1 q « 20s.However, when considering M 0 , m s 0 , R f 1 , R s 1 , and R x as unknown, the CRB increase to 192, 160s and 117, 183s, respectively.This difference can be understood with the geometric interpretation of the CRB that was introduced by Scharf and McWhorter (1993): The CRB of a single unknown parameter is simply the inverse squared ℓ 2 -norm of the signal's derivative wrt.respective parameter.In the case of multiple unknowns, the CRB is given by the inverse squared ℓ 2 -norm of the signal's orthogonalized derivative, i.e. after removing all components that are parallel to another gradient.Consequently, the key to a low CRB and, ultimately, a stable fit is to disentangle the individual derivatives such that their orthogonalized components are large.
The SIR sequence maps a bi-exponential recovery after a T 2 -selective inversion pulse.The relaxation curve is characterized by 5 parameters, which sets an upper limit for the number of model parameters.An additional spin disturbance increases the number of observations, which opens the door for the estimation of additional parameters.As illustrated above, however, the key to a stable estimation of additional parameters is to disentangle the signal's derivatives wrt. the unknown parameters.
The proposed hybrid-state approach resembles SIR in the T 2 -selective inversion pulse but adds a train of RF pulses for additional spin disturbances.The large number of RF pulses provides many degrees of freedom to optimize the spin trajectory for maximum disentanglement of all derivatives, as captured by the CRB.
The choice between models with different numbers of parameters generally entails a variance-bias trade-off.For the experimental design described in this paper, which was optimized for unconstrained qMT imaging, the CRB values associated with an unconstrained model, compared to a constraint model (R s 1 fixed to 2s), are increased only by a factor of 2 in R f 1 , which corresponds to an SNR decrease of ?2, assuming that the CRB is a tight bound.For all other parameters, the CRB increases by a factor smaller than 1.03, which highlights the effectiveness of the disentanglement.Comparing the CRB values to experimental designs that are optimized for constraint qMT (supporting Tab.S1 vs. Tab. 1 in Assländer et al. (2024)), the largest SNR penalty is a factor of 3 in R f 1 (CRB increase by a factor of 9) in comparison to SIR with a 4-parameter model (M 0 , m s 0 , R f 1 , and a B 1 proxy (Cronin et al., 2020)) and much less in all other parameters.These comparisons confirm the feasibility of unconstrained qMT with 9 parameters (plus a complex phase) with a manageable SNR penalty.

Constrained vs. unconstrained MT models
Our data confirms previous reports that estimates T s 1 !T f 1 for white matter at 3T (Helms and Hagberg, 2009;Gelderen et al., 2016;Manning et al., 2021;Samsonov and Field, 2021;Zaiss et al., 2022).We show that this finding has substantial implications for the estimation of the other model parameters.With a Taylor expansion, we show that T f 1 and m s 0 are underestimated if T s 1 " T f 1 is assumed (Section 2.1.1)and a comparison of our experimental data to the literature confirms this finding.Section 2.1.1 also highlights that the finding T s 1 !T f 1 implies that MT drives the observed longitudinal relaxation, not just immediately following RF irradiation but throughout the MR experiment: in such a spin system continuous magnetization transfer to the semi-solid spin pool is a key driver of the apparent T 1 relaxation.This stands in contrast to most MT literature, which assumes that the MT effect is mostly during RF irradiation and that, once the longitudinal magnetization of the two pools approaches each other (which happens at the time scale T x " 1{R x « 50ms), they relax independently.
Inserting unconstrained estimates of qMT parameters in white matter (Tab.2) into Eq.( 2) results in T f,a 1 " 1{R f,a 1 « 0.94s (supporting Tab.S2), which is consistent with mono-exponential estimates reported in the literature (T f,a 1 « 1.084s (Stanisz et al., 2005)).This concordance is expected for experiments with z f {m f 0 « z s {m s 0 , which can be achieved in inversion recovery experiments either by inverting both spin pools with a short RF-pulse (T RF !T s 2 )-which is not feasible in vivo, but was done by Stanisz et al. (2005) in their NMR experiments-or by choosing inversion times that fulfill T I " T x .Our pulse sequence does not fulfill this condition, which explains the deviating T f,a 1 « 1.429s when fitting a mono-exponential model to our data (Supporting Fig. S4).Koenig et al. (1990) suggested that myelin is the primary source of GM-WM contrast in T 1 -weighted MRI (Fig. 3c,d), an observation that extends to R f,a 1 maps (Fig. 4).In an unconstrained MT model, the pronounced GM-WM contrast shifts from R f 1 to m s 0 (Fig. 3e,f).This observation refines the finding of Koenig et al. (1990) by identifying MT as the primary mechanism that generates the observed GM-WM contrast.However, we also observe a subtle GM-WM contrast in T f 1 , which may suggest that myelin also facilitates direct longitudinal relaxation of the free spin pool beyond MT, possibly by interactions between water protons and the local magnetic field of myelin (or macromolecules in general, see Gossuin et al. (2000)).

Myelin as a contrast agent
This observation is consistent with our phantom experiments (Fig. 2), where R f 1 was found to be linearly dependent on the BSA concentration.

Iron as a contrast agent
R f 1 of the pallidum was shorter than that of all other ROIs analyzed in this study (Tab.2).Since iron is known to accumulate in the pallidum in the form of ferritin, this suggests a sensitivity of R f 1 to iron, which matches the reports by Vymazal et al. (1999) and Samsonov and Field (2021).Supporting Fig. S8 fits R f 1 as a function of the iron concentration in each ROI as taken from the literature, revealing a linear dependency (R 2 " 0.94).Repeating the same analysis for the transversal relaxation rate R f 2 reveals a much clearer linear dependency (R 2 " 0.9998), suggesting that T f 2 is more sensitive and specific to iron than R f 1 , in line with previous reports by Schenker et al. (1993); Vymazal et al. (1999); Gossuin et al. (2000).

Myelin water as a confounding factor
Our WM estimates of T f 2 deviate from previous reports (Stanisz et al., 2005).A possible explanation is that our model neglects contributions from myelin water (MW)or water trapped between the myelin sheaths-that has a characteristic T MW 2 « 10ms (Mackay et al., 1994).MW exchanges magnetization with myelin's macromolecular pool as well as the larger intra-/extra-axonal water pool, where the former exchange is faster than the latter (Stanisz et al., 1999;Manning et al., 2021).A saturation of the semi-solid pool could, thus, result in a saturation of the MW pool and, ultimately, its suppression.A subsequent estimate of the observed T f 2 -which comprises both the intra-/extraaxonal water pool and MW pool-would thus be dominated by the former and result in higher observed T f 2 values.By contrast, a CPMG sequence starts from thermal equilibrium and has more pronounced contributions from the MW pool, resulting in shorter observed T f 2 .However, a more detailed analysis is needed for a thorough understanding of these observed deviations.

Unconstrained qMT in multiple sclerosis
Supporting Fig. S5u highlights four MS lesions with a hypointense appearance in the MP-RAGE.Our data suggests that this hypointensity is primarily driven by a reduction of m s 0 , which was observed in most examined lesions (Fig. 7b).By contrast, we find that changes in the NAWM are primarily driven by T f 1 (Fig. 5).This juxtaposition of the different sources of contrast changes highlights the complexity of longitudinal relaxation in biological tissue.
In histology, MS lesions exhibit substantial heterogeneity in terms of varying degrees of remyelination, axonal damage, inflammation, and gliosis (Lassmann, 2018).Fig. 7 suggests that unconstrained qMT can delineate more independent information as compared to constrained qMT.Future work will aim to identify links between qMT parameters and pathological variability in MS lesions.
Another goal of this paper was to gauge the sensitivity of unconstrained qMT to subtle changes in normalappearing WM and GM that are not easily detectable with established (contrast-based) clinical sequences.We observed statistically significant deviations of T f 1 between individuals with MS and healthy controls, in particular, in the NAWM, which aligns with previous studies that performed mono-exponential T 1 -mapping (Vrenken et al., 2006a,c,b).Moreover, we found statistically significant deviations of T f 1 in subcortical GM structures.An analysis of NAWM in individuals with MS always bears the risk of contaminating the results with an incomplete exclusion of MS lesions or by voxels close to lesions.However, we have two reasons to believe that lesions and their surrounding tissue do not drive the observed changes in R f 1 .First, we predominantly observe changes in m s 0 in lesions, while m s 0 changes in NAWM are much less pronounced.Second, Vrenken et al. (2006b) demonstrated that the magnetization transfer ratio in NAWM changes with the distance to an MS lesion, but their mono-exponential T 1 estimates do not.Another limitation of this study is the small number of participants, which does not allow for adjustments, e.g., of the age difference between the two cohorts.Therefore, larger studies are needed to confirm this result.

Rapid, high-resolution qMT imaging
A major goal of this paper is to demonstrate the feasibility of unconstrained qMT imaging on a voxel-by-voxel basis.With a hybrid-state pulse sequence, we were able to extract unconstrained qMT maps with 1.24mm, 1.6mm, and 2.0mm isotropic resolution from 12min, 6min, and 4min scans, respectively.To the best of our knowledge, the presented maps are the first voxel-wise fits using an unconstrained MT model.However, we do observe a subtle blurring in our qMT maps compared to the MP-RAGE.The most likely cause is the smaller k-space coverage of the koosh-ball trajectory in comparison to a Cartesian trajectory: the koosh-ball trajectory with a nominal resolution of 1.0mm samples only the inner sphere of the 1.0mm kspace cube, similar to elliptical scanning, while the MP-RAGE samples the entire cube.We account for the reduced k-space coverage using the "effective" resolution of 1.24mm.Undersampling, regularized reconstruction, and incomplete motion correction might cause additional blurring.On the flip side, our image reconstruction models the spin dynamics, alleviating relaxation-induced blurring that is more prevalent in approaches like MP-RAGE (Mugler and Brookeman, 1990) or RARE (Hennig et al., 1986).

Future Directions
Our ongoing work includes clinical validation as well as efforts for further scan time reductions and improvements in resolution.To this end, we aim to replace the current RF pattern, which is a concatenation of separate optimizations, with a joint optimization of all unconstrained qMT parameters.Further, we are exploring more efficient kspace trajectories.Last, we anticipate that studies with the current pulse sequence will help identify the most clinically meaningful parameters.This information can then be fed back to our numerical optimization framework to optimize pulse sequences for more efficient estimation of these parameters.Optimizations of the sequence for particular parameters can be achieved using the employed CRB-based framework without imposing constraints on the parameters.

Conclusion
Our study builds on the work of Helms and Hagberg (2009) (Hallgren and Sourander, 1958).a: Except for WM and the thalamus, which were excluded from the fit, the longitudinal relaxation rate R f 1 " 1{T f 1 follows in good approximation a linear function (R 2 " 0.94).b: The transversal relaxation rate of the free pool R f 2 " 1{T f 2 is even better described by a linear model (R 2 " 0.9998). 1 " 1s, respectively, and the proposed unconstrained MT model.This analysis pools all normal-appearing white matter voxels of 4 healthy participants and 4 individuals with MS, respectively.Note that the T s 1 " T f 1 column depicts the same longitudinal relaxation time estimates twice to illustrate the differences to the unconstrained MT model.For all six WM and GM ROIs, we found statistically significant differences between participants with MS and controls.The markers ˚and ˚˚indicate statistically significant differences at the p ă 0.05 and p ă 0.01 levels.

Figure 2 :
Figure2: Phantom validation.Seven tubes filled with different concentrations of bovine serum albumin (BSA) were imaged at 1.5T and 2.9T.The box plots represent the median, 1 st and 3 rd quartile, and the whiskers the 1.5x the inter-quartile range or the maximum range, whichever is smaller.The median values of each tube's qMT estimates were fitted with a general linear model with the BSA concentration (c BSA ) and the field strength (B 0 ) as independent variables.The fitted coefficients are listed in Tab. 1.The brackets indicate outliers that were excluded from the GLM regression due to unstable parameter estimation of the semi-solid pool's characteristic, likely caused by the small pool size.

Fig. 4 Figure 3 :
Figure 3: Comparison of clinical contrasts (a-d) and quantitative magnetization transfer (qMT) maps (e-p) in a healthy volunteer.The qMT maps have an effective resolution of 1.24mm isotropic (acquired in 36min) compared to the 1mm isotropic of the clinical contrasts.We display here relaxation rates (R f,s 1,2 " 1{T f,s 1,2 ), where the superscripts f and s indicate the free and semi-solid pools, respectively.The size of the semi-solid spin pool is normalized by m s 0 `mf 0 " 1, and Rx denotes the exchange rate between the two pools.

Figure 4 :
Figure 4: Apparent quantitative MT maps when assuming T s 1 " T f 1 in a healthy volunteer.The maps were calculated voxel-wise with Eqs.(3)-(5) and based on the maps depicted in Fig. 3.Note the different color scale in R f,a 1 compared to Fig. 3.

Figure 5 :
Figure 5: Comparison of apparent qMT parameter estimates when assuming T s 1 " T f 1 (Eqs.(3)-(5)) to unconstrained parameter estimates.The box plots pool all normal-appearing white matter voxels of each participant.The markers ˚and ˚˚indicate statistically significant differences at the p ă 0.05 and p ă 0.01 levels in a comparison of each subject's median qMT parameter estimates between participants with MS and controls.

Figure 6 :
Figure6: ROI analysis of the unconstrained qMT model's parameters, pooled over all normal-appearing white matter (NAWM) voxels in each of the 4 individuals with MS and the 4 controls.The markers ˚and ˚˚indicate statistically significant differences at the p ă 0.05 and T s 1 in the anterior corpus callosum (p ă 0.01, p ă 0.05, p ă 0.05); • T f 1 in the posterior corpus callosum (p ă 0.01); • T f 1 in the cortical GM (p ă 0.05); • T f 1 and T f 2 in the caudate (p ă 0.05, p ă 0.05); • T f 1 in the pallidum (p ă 0.05); • T f 1 in the putamen (p ă 0.05).

Figure 7 :
Figure 7: Analysis of longitudinal relaxation in lesions pooled across all 4 participants with MS where each color corresponds to one individual.a The median size of the apparent semi-solid spin pool m s,a 0 vs the median apparent relaxation rate R f,a 1 .b Median m s 0 vs R f 1 as measured with the unconstrained MT model.The black arrows visualize the scaled eigenvectors of a PCA that quantify the independent variability in the respective model.

Figure 8 :
Figure 8: Quantitative MT maps of an individual with MS.The rows compare different (isotropic) effective resolutions that require different scan times.All scans were acquired with full brain coverage.The magnifications show a lesion to highlight the resolution differences and the consequent partial volume effects.

Figure S2 :
Figure S2: Parameter to noise ratio (PNR) of phantom scans at 2.9T.The PNR was measured voxel-by-voxel based on 25 repetitions of the same scan (reconstructed without any regularization) and the box plots show the inter-voxel variability of the estimated PNR.The x-markers identify the PNR as estimated by the median parameter estimates divided by the corresponding Cramèr-Rao bound (CRB).

Figure S3 :
Figure S3: Comparison of parameter estimates between unregularized and regularized reconstructions at 2.9T.Each box plot analyzes the parameter estimates in a sample, pooled over all voxels in respective ROI and over 25 repetitions of the scan.The utilized regularization, including the regularization parameter, was identical to the reconstruction of the in vivo data shown in Fig. 8.

Figure S6 :
Figure S6: Cramér-Rao bound (CRB) of the in vivo parameter estimates shown in Fig.S5a-f, normalized to resemble the inverse squared SNR of the parameter estimates (cf.Tab.S1).We used the unconstrained MT model and the corresponding estimates for the CRB calculations.The FLAIR image is provided for anatomical reference.

Figure S7 :
Figure S7: Apparent quantitative MT maps when assuming T s 1 " T f 1 in a participant with MS.The maps were calculated voxel-wise with Eqs.(3)-(5) and based on the maps depicted in Fig. S5.Note the different color scale in R f,a 1 compared to Fig. 3.

Figure S8 :
FigureS8: Relaxation rates, measured with our qMT approach, as a function of iron concentrations, taken from the literature(Hallgren and Sourander, 1958).a: Except for WM and the thalamus, which were excluded from the fit, the longitudinal relaxation rate R f 1 " 1{T f 1 follows in good approximation a linear function (R 2 " 0.94).b: The transversal relaxation rate of the free pool R f 2 " 1{T f 2 is even better described by a linear model (R 2 " 0.9998).

Figure S9 :
FigureS9: Comparison the parameter estimates between a Bloch model, two traditional MT models that assume T s 1 " T f 1 and T s 1 " 1s, respectively, and the proposed unconstrained MT model.This analysis pools all normal-appearing white matter voxels of 4 healthy participants and 4 individuals with MS, respectively.Note that the T s 1 " T f 1 column depicts the same longitudinal relaxation time estimates twice to illustrate the differences to the unconstrained MT model.

Figure S10 :
Figure S10: ROI analysis of the unconstrained model's T f1 .For all six WM and GM ROIs, we found statistically significant differences between participants with MS and controls.The markers ˚and ˚˚indicate statistically significant differences at the p ă 0.05 and p ă 0.01 levels.

Table 1 :
using the FLAIR images.The automated lesion segmentations were manually adjusted by FLR and ESB and subtracted from the ROI masks.After that, we eroded the outmost layer of voxels of each ROI a 0 a BSA a B0 p 1 {Tq a 2 p 1 {Tq Coefficients of a generalized linear model fit of the phantom data shown in Fig.