Unconstrained quantitative magnetization transfer imaging: disentangling $T_{1}$ of the free and semi-solid spin pools

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 Fig. 1, captures all protons bound in liquids where fast molecular motion causes an exponential relaxation of the transversal magnetization with a characteristic $T_{2}^{f}\gtrsim 50$ms (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 $\omega_{z}$, and the Rabi frequency of the RF pulses by $\omega_{y}$. For readability, we here use relaxation rates ($R_{1,2}^{f,s}=1/T_{1,2}^{f,s}$). 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_{2}^{s}\approx 10\upmu$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 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 $R_{2}^{s,l}(R_{2}^{s},\alpha,T_{\text{RF}})$. While exponential and non-exponential decays necessarily deviate, we can identify an $R_{2}^{s,l}$ that results in the same magnetization at the end of an RF pulse. To this end, $R_{2}^{s,l}$ depends on the flip angle $\alpha$ and the duration $T_{\text{RF}}$ of respective RF-pulse in addition to the biophysical parameter $R_{2}^{s}$. We neglect the $y^{s}$ component assuming, without loss of generality, $\omega_{x}=0$ and given that $R_{2}^{s,l}\gg\omega_{z}$. The exchange rate $R_{\text{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_{0}^{s}$ have been used. Here, we use $m_{0}^{f}+m_{0}^{s}=1$ so that $m_{0}^{s}$ 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 $\tilde{m}_{0}^{s}=m_{0}^{s}/m_{0}^{f}$ or, equivalently, normalize to $\tilde{m}_{0}^{f}=1$. The conversion between the two definitions is simply $m_{0}^{s}=\tilde{m}_{0}^{s}/(1+\tilde{m}_{0}^{s})$ and $\tilde{m}_{0}^{s}=m_{0}^{s}/(1-m_{0}^{s})$.

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_{\text{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 $\pi$ 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 $\pi$ between consecutive RF pulses. The pulses are separated by a $T_{\text{R}}=3.5\text{ms}$, which is approximately the minimal $T_{\text{R}}$ with which we can perform gradient encoding with $|k_{\max}|=\pi/1$mm 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 $\pi$ 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).

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 $m_{0}^{s}$, $R_{1}^{f}$, $R_{2}^{f}$, $R_{\text{x}}$, $R_{1}^{s}$, and $T_{2}^{s}$, while additionally accounting for $\omega_{z}$, $B_{1}^{+}=\omega_{y}/\omega_{y}^{\text{nominal}}$, 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_{0}^{s}=0.25$, $R_{1}^{f}=0.5/\text{s}$, $R_{2}^{f}=15.4/\text{s}$, $R_{\text{x}}=20/\text{s}$, $R_{1}^{s}=2$/s, $T_{2}^{s}=10\upmu\text{s}$, $\omega_{z}=0$, and $B_{1}^{+}=1$. The resulting spin trajectories (Fig. LABEL:fig:Spin_Dynamics) and the corresponding CRB values (Tab. LABEL:tab:CRB) are discussed in the supporting information Section LABEL:supsec:NumericalOptimizations. Supporting Section LABEL:supSec:Noise and Figs. LABEL:fig:Phantom_PNRCG, LABEL:fig:Phantom_repLLRvCG, and LABEL:fig:InVivo_CRB connect the CRB to experimental noise measurements.

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}|=\pi/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 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).

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\pm 8.7$; 3 female) and 4 healthy controls (age $28.8\pm 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.

  1.24mm isotropic in 12min

• 2.

  1.6mm isotropic in 6min

• 3.

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

We reconstructed the images with sub-space modeling (Liang, 2007; Huang et al., 2012; Christodoulou et al., 2018; McGivney et al., 2014; Tamir et al., 2017; Assländer et al., 2018; Zhao et al., 2018), i.e., we reconstructed coefficient images in the sub-space spanned by singular vectors from a coarse dictionary of signals (or fingerprints) (McGivney et al., 2014; Tamir et al., 2017; Assländer et al., 2018) and their orthogonalized gradients (Mao et al., 2023b). We used the optISTA algorithm (Jang et al., 2023), incorporating sensitivity encoding (Sodickson and Manning, 1997; Pruessmann et al., 2001) and locally low-rank regularization (Lustig et al., 2007; Trzasko and 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. B). 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; Nataraj et al., 2018; Duchemin et al., 2020; Zhang et al., 2022; Mao et al., 2023a). This approach includes a data-driven $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).

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 ($c_{\text{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_{0}^{s}=a_{\text{BSA}}\cdot c_{\text{BSA}}$, i.e., we can reject neither the null hypothesis that $m_{0}^{s}$ is independent of $B_{0}$ nor that it vanishes at $c_{\text{BSA}}=0$.

For $R_{1}^{f}$, our data supports the linear model $R_{1}^{f}=a_{0}+a_{\text{BSA}}\cdot c_{\text{BSA}}+a_{2}\cdot c_{\text{BSA}}\cdot B_{0}$. I.e., the data suggests a linear dependence of $R_{1}^{f}$ on the BSA concentration where the slope also depends on the field strength. By contrast, we did not observe a dependency of $R_{1}^{f}$ of pure water ($c_{\text{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 $\tau_{c}=5\text{ps}$; (Bloembergen et al., 1948)). Further, the intercept $a_{0}=0.28/\text{s}$ with a 95% confidence interval $[0.15,0.41]/\text{s}$ agrees with the $0.255/\text{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_{2}^{f}$, our data is consistent with a linear dependence on the BSA concentration and no dependence on $B_{0}$ ($R_{2}^{f}=a_{\text{BSA}}\cdot c_{\text{BSA}}$). 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/\text{s}$ is not physical, the 95% confidence interval $[-7.6,1.8]/\text{s}$ includes the $R_{2}^{f}\approx 0.255/\text{s}$ predicted by the BPP theory.

We detected no dependence of $R_{1}^{s}$ on $c_{\text{BSA}}$, but a statistically significant dependence on $B_{0}$ ($R_{1}^{s}=a_{0}+a_{B_{0}}\cdot B_{0}$), which is consistent with the reports of Wang et al. (2020). For $R_{\text{x}}$, we observe neither a $c_{\text{BSA}}$ nor a $B_{0}$ dependency, and for $T_{2}^{s}$, the $B_{0}$ dependency is just above the 95% significance threshold.

Beyond the linear model, we observe increased variability of $R_{\text{x}}$, $R_{1}^{s}$, and $T_{2}^{s}$ estimates with decreasing $c_{\text{BSA}}$, which likely stems from the smaller spin pool size. For this reason, we excluded the most extreme case ($c_{\text{BSA}}=0.05$) from the GLM fits as indicated by the brackets in Fig. 2.

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_{0}^{s}$, $R_{1}^{f}$, and $R_{2}^{f}$. 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_{\text{x}}$ and $R_{1}^{s}$ maps exhibit reduced image quality, consistent with their higher CRB values (Tab. LABEL:tab:CRB). Also consistent with its large CRB, the $T_{2}^{s}$ 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_{2}^{f}$ (Fig. 3i,j, and Fig. LABEL:fig:InVivo_MSc) 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_{0}^{s}$, followed by $R_{1}^{f}$. In $R_{2}^{f}$, however, we observe only a subtle contrast between cortical GM and WM. An ROI analysis confirms this finding: we estimated $T_{2}^{f}=(83\pm 15)$ms and $(76.9\pm 8.3)$ms for cortical GM and WM, respectively, which is a smaller difference compared to the difference between previously reported values ($(99\pm 7)$ vs. $(69\pm 3)$ms; (Stanisz et al., 2005)). Consistent with previous reports, we observe the shortest $T_{2}^{f}=(59.3\pm 5.6)$ms in the pallidum (Fig. 3i). The exchange rate $R_{\text{x}}$, $R_{1}^{s}$, and $T_{2}^{s}$ exhibit little GM-WM contrast and we note that the most prominent contrast in $R_{\text{x}}$ and $R_{1}^{s}$ 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_{0}^{s}$ 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.

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.

Many pulse sequences are sensitive to $R_{1}^{f}$ and $R_{1}^{s}$. Using SIR (Dortch et al., 2018) as an example, this can be illustrated with normalized CRB values (cf. supporting Tab. LABEL:tab:CRB) under the assumption that only the respective parameter is unknown: norm. $\text{CRB}(R_{1}^{f})\approx 25\text{s}$ and norm. $\text{CRB}(R_{1}^{s})\approx 20\text{s}$. However, when considering $M_{0}$, $m_{0}^{s}$, $R_{1}^{f}$, $R_{1}^{s}$, and $R_{x}$ as unknown, the CRB increase to $192,160\text{s}$ and $117,183\text{s}$, 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 $\ell_{2}$-norm of the signal’s derivative wrt. respective parameter. In the case of multiple unknowns, the CRB is given by the inverse squared $\ell_{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_{1}^{s}$ fixed to 2s), are increased only by a factor of 2 in $R_{1}^{f}$, which corresponds to an SNR decrease of $\sqrt{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. LABEL:tab:CRB vs. Tab. 1 in Assländer et al. (2024)), the largest SNR penalty is a factor of 3 in $R_{1}^{f}$ (CRB increase by a factor of 9) in comparison to SIR with a 4-parameter model ($M_{0}$, $m_{0}^{s}$, $R_{1}^{f}$, 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.

Our data confirms previous reports that estimates $T_{1}^{s}\ll T_{1}^{f}$ 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_{1}^{f}$ and $m_{0}^{s}$ are underestimated if $T_{1}^{s}=T_{1}^{f}$ 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_{1}^{s}\ll T_{1}^{f}$ 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_{\text{x}}=1/R_{\text{x}}\approx 50$ms), they relax independently.

Inserting unconstrained estimates of qMT parameters in white matter (Tab. 2) into Eq. (2) results in $T_{1}^{f,a}=1/R_{1}^{f,a}\approx 0.94$s (supporting Tab. LABEL:tab:ROImean_apparent), which is consistent with mono-exponential estimates reported in the literature ($T_{1}^{f,a}\approx 1.084$s (Stanisz et al., 2005)). This concordance is expected for experiments with $z^{f}/m_{0}^{f}\approx z^{s}/m_{0}^{s}$, which can be achieved in inversion recovery experiments either by inverting both spin pools with a short RF-pulse ($T_{\text{RF}}\ll T_{2}^{s}$)—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_{\text{I}}\gg T_{\text{x}}$. Our pulse sequence does not fulfill this condition, which explains the deviating $T_{1}^{f,a}\approx 1.429$s when fitting a mono-exponential model to our data (Supporting Fig. LABEL:fig:InVivo_cf_Models).

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_{1}^{f,a}$ maps (Fig. 4). In an unconstrained MT model, the pronounced GM-WM contrast shifts from $R_{1}^{f}$ to $m_{0}^{s}$ (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_{1}^{f}$, 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 $R_{1}^{f}$ was found to be linearly dependent on the BSA concentration.

$R_{1}^{f}$ 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_{1}^{f}$ to iron, which matches the reports by Vymazal et al. (1999) and Samsonov and Field (2021). Supporting Fig. LABEL:fig:Iron fits $R_{1}^{f}$ 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_{2}^{f}$ reveals a much clearer linear dependency ($R^{2}=0.9998$), suggesting that $T_{2}^{f}$ is more sensitive and specific to iron than $R_{1}^{f}$, in line with previous reports by Schenker et al. (1993); Vymazal et al. (1999); Gossuin et al. (2000).

Our WM estimates of $T_{2}^{f}$ 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_{2}^{\text{MW}}\approx 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 (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_{2}^{f}$—which comprises both the intra-/extra-axonal water pool and MW pool—would thus be dominated by the former and result in higher observed $T_{2}^{f}$ values. By contrast, a CPMG sequence starts from thermal equilibrium and has more pronounced contributions from the MW pool, resulting in shorter observed $T_{2}^{f}$. However, a more detailed analysis is needed for a thorough understanding of these observed deviations.

Supporting Fig. LABEL:fig:InVivo_MSu 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_{0}^{s}$, which was observed in most examined lesions (Fig. 7b). By contrast, we find that changes in the NAWM are primarily driven by $T_{1}^{f}$ (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 normal-appearing WM and GM that are not easily detectable with established (contrast-based) clinical sequences. We observed statistically significant deviations of $T_{1}^{f}$ 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_{1}^{f}$ 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_{1}^{f}$. First, we predominantly observe changes in $m_{0}^{s}$ in lesions, while $m_{0}^{s}$ 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.

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