Unconstrained quantitative magnetization transfer imaging: Disentangling T₁ of the free and semi-solid spin pools Open Access

Longitudinal relaxation is a vital contrast mechanism in magnetic resonance imaging (MRI). For example, the MP-RAGE (Mugler & Brookeman, 1990) pulse sequence generates excellent gray matter (GM)–white matter (WM) contrast and—compared to mostly $T2$-weighted pulse sequences like FLAIR (Hajnal et al., 1992)—may be more specific to the underlying tissue changes in multiple sclerosis (MS) lesions (Bagnato et al., 2003; Barkhof, 1999). 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 & 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 (⁠$T2s≈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 different parameters. In order to mitigate the consequent noise amplification, most quantitative MT (qMT) approaches constrain $T1s=T1f$ resulting in $T1s≈1.1$s (Dortch et al., 2011; Yarnykh, 2002) or, similarly, simply assume that $T1s=1$s (Henkelman et al., 1993; Morrison & Henkelman, 1995). However, recent studies have suggested that $T1s≈0.3$s and $T1f≈2$s in white matter at 3T (Gelderen et al., 2016; Helms & Hagberg, 2009; Manning et al., 2021; Samsonov & Field, 2021). These studies overcame the challenges of an unconstrained model by either using brain-wide estimates of $T1s$ and/or $T1f$ (Gelderen et al., 2016; Samsonov & Field, 2021) or fitting the MT model to NMR samples (Manning et al., 2021) or a single large ROI averaged over multiple participants (Helms & 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 encoding 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 4min, respectively.

We use the MT model described in Assländer et al. (2022, 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 Figure 1, captures all protons bound in liquids where fast molecular motion causes an exponential relaxation of the transversal magnetization with a characteristic $T2f≳50$ms (Bloembergen et al., 1948). The free pool’s magnetization is described by the Cartesian coordinates $xf$⁠, $yf$⁠, $zf$⁠, 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 (⁠$R1,2f,s=1​/T1,2f,s$⁠). The magnetization components $xs,zs$ of the semi-solid spin pool, sketched in purple in Figure 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 $T2s≈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 & Henkelman, 1995). The non-exponential 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 $R2s,l(R2s,α,TRF)$⁠. While exponential and non-exponential decays necessarily deviate, we can identify an $R2s,l$ that results in the same magnetization at the end of an RF pulse. To this end, $R2s,l$ depends on the flip angle $α$ and the duration $TRF$ of respective RF-pulse in addition to the biophysical parameter $R2s$⁠. We neglect the $ys$ component assuming, without loss of generality, $ωx=0$ and given that $R2s,l≫ωz$⁠. The exchange rate $Rx$ 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 $m0s$ have been used. Here, we use $m0f+m0s=1$ so that $m0s$ describes the fraction of the overall spin pool, a definition which has also been used, for example, by Yarnykh (2012). Other papers, such as Henkelman et al. (1993); Gochberg and Gore (2003), measure $m˜0s=m0s​/​m0f$ or, equivalently, normalize to $m˜0f=1$⁠. The conversion between the two definitions is simply $m0s=m˜0s​/(1+m˜0s)$ and $m˜0s=m0s​/(1−m0s)$⁠.

In the absence of RF pulses (⁠$ωy=0$⁠), we can isolate the longitudinal components of Eq. (1):

An eigendecomposition of the Hamiltonian in Eq. (2) has three distinct eigenvalues (Gochberg & Gore, 2003; Henkelman et al., 1993; 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 $R1f,a$ that is approximated by the following Taylor expansion at $R1s=R1f$⁠:

The MT contributions to $R1f,a$ therefore depend foremost on the macromolecular pool size $m0s$ and the two relaxation rates. Higher order terms additionally depend on the exchange rate $Rx$⁠. Eq. (3) shows that $R1f,a≈R1f$ is a reasonable approximation only if $m0s(R1s−R1f)≪R1f$⁠. Otherwise, this linear correction term contributes significantly to $R1f,a$⁠, making MT an important driver of longitudinal relaxation. For example, let us assume $R1f=0.5​/$s, $R1s=3$/s, and $m0s=0.2$⁠. In this case, the linear correction term is 0.5/s and, thus, $R1f,a≈1.0​/​s≈⁄R1f$⁠.

The largest eigenvalue (in absolute value) is given by

which is dominated by the exchange rate $Rx$ 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):

Eq. (5) reveals that $m0s$ is underestimated when assuming $R1s=R1f$⁠. To give a sense of the magnitude of this bias, we can insert the above example values and further assume $Rx=15​/$s, which results in $m0s,a≈0.15$ instead of the underlying $m0s=0.2$⁠.

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 $TR$⁠. 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 & 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 $TR=3.5ms$⁠, which is approximately the minimal $TR$ with which we can perform gradient encoding with $|kmax|=π​/1$mm and avoid stimulating the peripheral nerves. After 1142 RF-pulses, that is, 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 $T2$-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, 2019b), 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).

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) (Cramér, 1946, Rao, 1945), 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 & McFadden, 1994; Wu, 1981). Nonetheless, the CRB can be used as a proxy for the “SNR-efficiency” or “conditioning” (Haldar & Kim, 2019; Zhao et al., 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 $m0s$⁠, $R1f$⁠, $R2f$⁠, $Rx$⁠, $R1s$⁠, and $T2s$⁠, while additionally accounting for $ωz$⁠, $B1+=ωy​/​ωynominal$⁠, and the scaling factor $M0$ as unknowns, where $M0$ 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 $m0s=0.25$⁠, $R1f=0.5​/​s$⁠, $R2f=15.4​/​s$⁠, $Rx=20​/​s$⁠, $R1s=2$/s, $T2s=10μs$⁠, $ωz=0$⁠, and $B1+=1$⁠. The resulting spin trajectories (Supporting Figure S1) and the corresponding CRB values (Supporting Table S1) are discussed in the Supporting Information Section S1. Supporting Section S2 and Figures S2, S3, and S6 connect the CRB to experimental noise measurements.

We built a custom phantom composed of cylindrical 50 mL 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 10min.

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 $|kmax|=π​/1$mm). 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 regularized reconstruction. We changed the direction of the k-space spokes with a 2D golden means pattern (Chan et al., 2009; Winkelmann et al., 2007) that is reshuffled to improve the k-space coverage for each time point and to minimize eddy current artifacts (Flassbeck & Assländer, 2023).

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 relapsing-remitting 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:

• 1.24mm isotropic in 12min

• 1.6mm isotropic in 6min

• 2.0mm isotropic in 4min.

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 basic 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.

We reconstructed the images with sub-space modeling (Christodoulou et al., 2018; Huang et al., 2012; Liang, 2007; McGivney et al., 2014; Tamir et al., 2017; Zhao et al., 2018), that is, we reconstructed coefficient images in the sub-space spanned by singular vectors from a coarse dictionary of signals (or fingerprints) (Assländer et al., 2018; McGivney et al., 2014; Tamir et al., 2017) and their orthogonalized gradients (Mao et al., 2023). We used the optISTA algorithm (Jang et al., 2023), incorporating sensitivity encoding (Pruessmann et al., 2001; Sodickson & Manning, 1997) and locally low-rank regularization (Lustig et al., 2007; Trzasko & Manduca, 2011; Zhang et al., 2015) to reduce residual undersampling artifacts and noise. We implemented this reconstruction in Julia and made the source code publicly available (cf. Section Data and Code Availability). A more detailed description of the reconstruction can be found in Tamir et al. (2017) and Assländer et al. (2018).

For the 36min in vivo scans, we used separate sub-spaces for each 6min sub-scan, implemented as a block-diagonal 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.

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; Duchemin et al., 2020; Mao et al., 2024; Nataraj et al., 2018; Zhang et al., 2022). This approach includes a data-driven $B0$ and $B1+$ 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 Figure 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. (2024).

For the 36min reference scans, we registered the skull-stripped (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, 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) 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 to reduce partial volume effects with other tissues and to ensure that all ROI voxels are at least one voxel away from any lesion.

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 (⁠$cBSA$⁠) and the magnetic field strength (⁠$B0$⁠) as independent variables (Fig. 2; Table 1).

The estimates of the semi-solid spin pool size are consistent with the linear model $m0s=aBSA·cBSA$⁠, that is, we can reject neither the null hypothesis that $m0s$ is independent of $B0$ nor that it vanishes at $cBSA=0$⁠.

For $R1f$⁠, our data support the linear model $R1f=a0+aBSA·cBSA+a2·cBSA·B0$⁠, that is, the data suggest a linear dependence of $R1f$ on the BSA concentration where the slope also depends on the field strength. By contrast, we did not observe a dependency of $R1f$ of pure water (⁠$cBSA=0$⁠) on $B0$⁠. 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 $a0=0.28​/​s$ with a 95% confidence interval $[0.15,0.41]​/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 $R2f$⁠, our data are consistent with a linear dependence on the BSA concentration and no dependence on $B0$ (⁠$R2f=aBSA·cBSA$⁠). The latter is also consistent with the very small change predicted by the BPP theory (approximately 0.001%). While the negative intercept $a0=−2.9​/​s$ is not physical, the 95% confidence interval $[−7.6,1.8]​/​s$ includes the $R2f≈0.255​/​s$ predicted by the BPP theory.

We detected no dependence of $R1s$ on $cBSA$⁠, but a statistically significant dependence on $B0$ (⁠$R1s=a0+aB0·B0$⁠), which is consistent with the reports of Wang et al. (2020). For $Rx$⁠, we observe neither a $cBSA$ nor a $B0$ dependency, and for $T2s$⁠, the $B0$ dependency is just above the 95% significance threshold.

Beyond the linear model, we observe increased variability of $Rx$⁠, $R1s$⁠, and $T2s$ estimates with decreasing $cBSA$⁠, which likely stems from the smaller spin pool size. For this reason, we excluded the most extreme case (⁠$cBSA=0.05$⁠) from the GLM fits as indicated by the brackets in Figure 2.

Figure 3 demonstrates the feasibility of unconstrained qMT imaging with a hybrid-state pulse sequence, that is, 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 $m0s$⁠, $R1f$⁠, and $R2f$⁠. 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 $Rx$ and $R1s$ maps exhibit reduced image quality, consistent with their higher CRB values (Supporting Table S1). Also consistent with its large CRB, the $T2s$ map has the highest noise and artifact levels, which might also be, in part, due to a residual $B1+$ artifact caused by incomplete spoiling of the inversion pulse. We also find subtle residual $B1+$ artifacts in $R2f$ (Fig. 3i, j; Supporting Fig. S5c) and residual $B0$ 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 $B0$ and $B1+$ correction.

Among all qMT parameters, we observe the largest quantitative GM-WM contrast in $m0s$⁠, followed by $R1f$⁠. In $R2f$⁠, however, we observe only a subtle contrast between cortical GM and WM. An ROI analysis confirms this finding: we estimated $T2f=(83±15)$ms and $(76.9±8.3)$ms for cortical GM and WM, respectively, which is a smaller difference compared to the difference between previously reported values (⁠$(99±7)$ms vs. $(69±3)$ms; (Stanisz et al., 2005)). Consistent with previous reports, we observe the shortest $T2f=(59.3±5.6)$ms in the pallidum (Fig. 3i). The exchange rate $Rx$⁠, $R1s$⁠, and $T2s$ exhibit little GM-WM contrast and we note that the most prominent contrast in $Rx$ and $R1s$ occurs in voxels subject to partial volume effects and in CSF. In voxels with partial volume, the model might be inaccurate and the small $m0s$ 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 Table 2 for selected WM and GM structures.

Figure 4 depicts the apparent qMT parameters associated with a $T1s=T1f$ constrained model (Eqs. (3)–(5)). Figure 5 compares the apparent qMT parameters of those fitted with the unconstrained model, with more brain ROIs analyzed in Supporting Table S2. Below, we discuss the salient differences in white and cortical gray matter.

With the unconstrained model, we estimated a substantially different $T1f=(1.84±0.17)$s from $T1s=(0.34±0.10)$s. Using Eq. (3) to calculate the apparent $T1f,a$⁠, we estimate $T1f,a=(0.941±0.069)s$⁠, which approximately matches literature values ((⁠$1.084±0.045$⁠)s (Stanisz et al., 2005)).

With the unconstrained MT model, we estimated $m0s=0.212±0.022$⁠, consistent with literature estimates using the same model (⁠$0.172±0.043$⁠; (Helms & Hagberg, 2009)). The apparent pool size (Eq. (5)) is $m0s,a=0.151±0.022$⁠, which also matches the constrained estimates in the literature (⁠$0.139±0.028$⁠; (Stanisz et al., 2005) and $0.118±0.050$⁠; (Helms & Hagberg, 2009)).

The exchange rate estimated with the unconstrained MT model, $Rx=(13.6±1.1)$/s, is slightly lower compared to the corresponding literature ((⁠$18.1±3.6$⁠)/s (Helms & Hagberg, 2009)) as is the the apparent exchange rate (Eq. (4)): $Rxa=(16.1±1.2)$/s compared to (⁠$23±4$⁠)/s (Stanisz et al., 2005). Notwithstanding, our analysis matches previous findings that the constraint $T1s=T1f$ biases $Rx$ to larger values (Helms & Hagberg, 2009).

Supporting Figures 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.

An ROI analysis of the cortical gray matter, averaged over all healthy volunteers, reveals trends similar to the WM analysis: $T1f=(2.46±0.56)$s and $T1s=(0.42±0.40)$s differ substantially from one another. The apparent $T1f,a=(1.62±0.23)s$ is in line with the mono-exponential estimate (⁠$1.82±0.11$⁠)s measured by Stanisz et al. (2005).

As expected, the estimated $m0s=0.098±0.026$ is both smaller than that measured in WM and similar to the literature value $m0s=0.086$ (derived from $m˜0s=0.094$⁠) estimated with the unconstrained MT model (Helms et al., 2004). The estimated $m0s,a=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 $Rx=(14.0±3.1)$/s as well as $Rxa=(16.4±3.4)$/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)).

Figure 6 compares all 6 unconstrained qMT parameters in an ROI spanning the entire NAWM between individuals with RRMS and healthy controls. We observe the most distinct differences in $T1f$⁠: the median $T1f$ across the NAWM of each MS subject averaged over all participants with MS was 98ms larger than in controls (⁠$p<0.01$⁠). By comparison, the apparent $T1f,a$ (Eq. (3)) differs only by 19ms (⁠$p<0.05$⁠; cf. Fig. 5).

In NAWM, the median $T2f$ 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 Table 2, we found statistically significant changes of

• $T1f$⁠, $T2f$⁠, and $T1s$ in the anterior corpus callosum (⁠$p<0.01$⁠, $p<0.05$⁠, and $p<0.05$⁠);

• $T1f$ in the posterior corpus callosum (⁠$p<0.01$⁠);

• $T1f$ in the cortical GM (⁠$p<0.05$⁠);

• $T1f$ and $T2f$ in the caudate (⁠$p<0.05$ and $p<0.05$⁠);

• $T1f$ in the pallidum (⁠$p<0.05$⁠);

• $T1f$ in the putamen (⁠$p<0.05$⁠).

With the constrained MT model, we only found significant differences in $T1f,a$ of the putamen (⁠$p<0.05$⁠). Supporting Figure S10 depicts $T1f$ for the ROIs above.

In MS lesions, we observe a substantial reduction of $m0s$ (Fig. 8; Supporting Fig. S5) and $m0s,a$ (Supporting Fig. S7) relative to the NAWM, consistent with the expected demyelination.

When jointly analyzing $T1f,a$ and $m0s,a$ 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.

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

Many pulse sequences are sensitive to $R1f$ and $R1s$⁠. Using SIR (Dortch et al., 2018) as an example, this can be illustrated with normalized CRB values (cf. Supporting Table S1) under the assumption that only the respective parameter is unknown: norm. $CRB(R1f)≈25s$ and norm. $CRB(R1s)≈20s$⁠. However, when considering $M0$⁠, $m0s$⁠, $R1f$⁠, $R1s$⁠, and $Rx$ 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, that is, 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 $T2$-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 $T2$-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 constrained model (⁠$R1s$ fixed to 2s), are increased only by a factor of 2 in $R1f$⁠, 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 constrained qMT (Supporting Table S1 vs. Table 1 in Assländer et al. (2024)), the largest SNR penalty is a factor of 3 in $R1f$ (CRB increase by a factor of 9) in comparison to SIR with a 4-parameter model (⁠$M0$⁠, $m0s$⁠, $R1f$⁠, and a $B1+$ 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.

Our data confirm previous reports that estimates $T1s≪T1f$ for white matter at 3T (Gelderen et al., 2016; Helms & Hagberg, 2009; Manning et al., 2021; Samsonov & 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 $T1f$ and $m0s$ are underestimated if $T1s=T1f$ is assumed (Section 2.2) and a comparison of our experimental data to the literature confirms this finding. Section 2.2 also highlights that the finding $T1s≪T1f$ 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 $T1$ relaxation. This stands in contrast to most MT literature, which assumes that MT effects the free spin pool mostly during RF irradiation and that, once the longitudinal magnetization of the two pools approaches each other (which happens at the time scale $Tx=1​/​Rx≈50$ms), the two pools relax independently.

Inserting unconstrained estimates of qMT parameters in white matter (Table 2) into Eq. (2) results in $T1f,a=1​/R1f,a≈0.94$s (Supporting Table S2), which is consistent with mono-exponential estimates reported in the literature (⁠$T1f,a≈1.084$s (Stanisz et al., 2005)). This concordance is expected for experiments with $zf​/​m0f≈zs​/​m0s$⁠, which can be achieved in inversion recovery experiments either by inverting both spin pools with a short RF-pulse (⁠$TRF≪T2s$⁠)—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 $TI≫Tx$⁠. Our pulse sequence does not fulfill this condition, which explains the deviating $T1f,a≈1.429$s 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 $T1$-weighted MRI (Fig. 3c, d), an observation that extends to $R1f,a$ maps (Fig. 4b). In an unconstrained MT model, the pronounced GM-WM contrast shifts from $R1f$ to $m0s$ (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 $T1f$⁠, 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)). This observation is consistent with our phantom experiments (Fig. 2), where $R1f$ was found to be linearly dependent on the BSA concentration.

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

Our WM estimates of $T2f$ 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 $T2MW≈10$ms (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 (Manning et al., 2021; Stanisz et al., 1999). 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 $T2f$—which comprises both the intra-/extra-axonal water pool and MW pool—would thus be dominated by the former and result in higher observed $T2f$ values. By contrast, a CPMG sequence starts from thermal equilibrium and has more pronounced contributions from the MW pool, resulting in shorter observed $T2f$⁠. However, a more detailed analysis is needed for a thorough understanding of these observed deviations.

Supporting Figure S5u highlights four MS lesions with a hypointense appearance in the MP-RAGE. Our data suggest that this hypointensity is primarily driven by a reduction of $m0s$⁠, which was observed in most examined lesions (Fig. 7b). By contrast, we find that changes in the NAWM are primarily driven by $T1f$ (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). Figure 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 normal-appearing WM and GM that are not easily detectable with established (contrast-based) clinical sequences. We observed statistically significant deviations of $T1f$ between individuals with MS and healthy controls, in particular, in the NAWM, which aligns with previous studies that performed mono-exponential $T1$-mapping (Vrenken et al., 2006a, 2006b, 2006c). Moreover, we found statistically significant deviations of $T1f$ 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 $R1f$⁠. First, we predominantly observe changes in $m0s$ in lesions, while $m0s$ 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 $T1$ estimates do not. Another limitation of this study is the small number of participants, which does not allow for adjustments, for example, of the age difference between the two cohorts. Therefore, larger studies are needed to confirm this result.

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 k-space 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 & Brookeman, 1990) or RARE (Hennig et al., 1986).

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 k-space 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.

Our study builds on the work of Helms and Hagberg (2009); Gelderen et al. (2016); Manning et al. (2021); Samsonov and Field (2021); and Zaiss et al. (2022), who pioneered unconstrained fitting with Henkelman’s two-pool magnetization transfer model. By utilizing the encoding power of the hybrid state (Assländer et al., 2019b), we improved the sensitivity of the MRI data to the model’s parameters, enabling an unconstrained fit of the MT model to each voxel separately. Our results confirm previous observations of the substantially different longitudinal relaxation times of the free and semi-solid spin pools. The results also suggest a potential clinical value of unconstrained qMT for individuals with MS via the detection of changes in the NAWM and the characterization of MS lesions.

In order to promote reproducibility, we provide the latest version (v0.8.0, DOI: 10.5281/zenodo.7433494) of the sequence optimization and signal simulation source code on https://github.com/JakobAsslaender/MRIgeneralizedBloch.jl. They are written in the open-source language Julia and registered to the package manager as “MRIgeneralizedBloch.jl.” The package documentation and tutorials can be found at https://JakobAsslaender.github.io/MRIgeneralizedBloch.jl. The tutorials render the code in HTML format with interactive figures and links to Jupyter notebooks that can be launched in a browser without local installations using binder. We also provide the source code of the image reconstruction package at https://github.com/JakobAsslaender/MRFingerprintingRecon.jl. For the presented data, we used v0.6.0. The qMT maps of all participants are available at https://zenodo.org/records/10729741.

Jakob Assländer: Conceptualization; Data curation; Formal analysis; Funding acquisition; Investigation; Methodology; Project administration; Software; Supervision; Visualization; Writing—original draft; and Writing—review & editing. Andrew Mao: Funding acquisition; Software; and Writing—review & editing. Elisa Marchetto: Motion Correction; Writing—review & editing. Erin S. Beck: Methodology (lesion segmentation); Conceptualization; and Writing—review & editing. Francesco La Rosa: Methodology (lesion segmentation); Writing—review & editing. Robert W. Charlson: Resources (patient recruitment); Writing—review & editing. Timothy Shepherd: Conceptualization; Writing—review & editing. Sebastian Flassbeck: Data curation; Investigation; Software; and Writing—review & editing.

This research was supported by the NIH/NINDS grant R01 NS131948, NIH/NIBIB grant R21 EB027241, and was performed under the rubric of the Center for Advanced Imaging Innovation and Research, an NIBIB Biomedical Technology Resource Center (NIH P41 EB017183). AM acknowledges support from the NIH/NIA (T32 GM136573 and F30 AG077794). FLR is supported by the Swiss National Science Foundation (SNF) Postdoc Mobility Fellowship (P500PB_206833).

None of the authors have any competing interests to declare.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00177.