We implemented a repeated IR spin-echo EPI pulse sequence with IDE preparation that allows independent control of T1 and MD weightings (Figure 1). The sequence is used to acquire whole-brain DWIs with a conventional single-shot EPI readout trajectory to prevent potential imaging artifacts arising from motion-induced phase inconsistencies in multi-shot in vivo acquisitions of diffusion-prepared signals.

The IR-preparation consists of a slice-selective adiabatic radio-frequency (RF) inversion pulse with a hyperbolic secant amplitude envelope and a 10.2 ms duration, followed by a gradient crusher. The IR-preparation and corresponding slice acquisition modules are interleaved to maintain the same TI and TR values for all slices (Park et al., 1985; Oh et al., 1991; Listerud et al., 1996). For any given slice, the slice-selective IR RF pulse excites the same location as the corresponding excitation and refocusing RF pulses but uses a 1.25 larger slice thickness to mitigate the loss of inversion efficiency due to mismatch of slice profiles. Moreover, to prevent cross-talk between IR pulses in adjacent slices, even and odd slices can be acquired in separate packets, with interleaved slice acquisition orders within each packet. The duration between two consecutive IR modules (and slice acquisitions) is determined by the TR and the number of slices per packet. In each scan, the T1-weighting is determined by fixing both TI and TR, while in consecutive TRs the diffusion-weighting is adjusted to acquire IR-DWIs with different b-values. The complete experiment consists of several repeated IR-prepared EPI scans, each with a different set of fixed (TI, TR) parameters determining different T1-weightings. To account for T1-encoding, both TI and TR are explicitly incorporated into the signal equation, as described in the following section.

The IDE-preparation module consists of a numerically optimized diffusion-weighting gradient waveform (Sjölund et al., 2015; Avram et al., 2019) that provides isotropic diffusion encoding in a single excitation (Mori and Van Zijl, 1995; Wong et al., 1995), but with a more efficient allocation of clinical diffusion gradient pulses before and after the 180° RF pulse (Avram et al., 2019). The desired maximum b-value determines the timings of the diffusion gradient pulses, and consequently the minimum TE. Rotation-invariant measurements with different b-values are acquired by simply scaling down the diffusion gradient pulse amplitudes. If we assume that the voxel is composed of multiple water pools, each with its own microscopic diffusion tensor whose corresponding diffusion ellipsoid has an arbitrary size, shape, and orientation, then the IDE preparation provides isotropic weighting by the Trace, or MD, for each of these microscopic tensors (Avram et al., 2019), allowing a spectroscopic quantitation of subvoxel MDs. Images with a fixed T1-weighting (TI and TR) but different b-values can be acquired efficiently in consecutive TRs of the same repeated IR scan. Because the T1 and diffusion weightings are encoded in the longitudinal and transverse magnetizations respectively, we can acquire signals with an arbitrary range of joint TI-b parameters. This contrasts with T2-diffusion CS-MRI measurements, where the practical range of joint TE-b parameters is limited by the minimum achievable TE as a function of b-value.

For a repeated IR spin-echo experiment (Figure 1) with perfect adiabatic inversion, the signal attenuation due to T1 and diffusion weightings in a homogeneous spin system with longitudinal relaxation time constant T₁ and mean diffusivity D̄ is:

By systematically varying the weighting parameters we can efficiently acquire multiple magnitude IR-IDE DWIs with a wide range of combined T1-MD weightings.

Subject and physiological motions can produce undesired phase contributions in DWIs acquired in vivo. Therefore, we can only reliably measure magnitude IR-IDE DWIs. However, in order to discriminate a superposition of longitudinal inversion recovery processes from subvoxel water pools with distinct T1s it is necessary to use phase-sensitive IR data that can capture both positive and negative values of longitudinal magnetization within a voxel. We adapted a simple method (Bakker et al., 1984) for estimating the polarity of IR measurements from multiple measurements with different TIs. Specifically, for a series of magnitude measurements |M(TI_k)|, associated with the increasing sequence of inversion times, TI_k, we first identify the minimum magnitude signal |M_m| = min|M(TI_k)| and its corresponding TI_m. Next, we invert the polarity of all magnitude measurements |M_k| with TI_k < TI_m and remove the minimum point at TI_m, to create a set of phase-sensitive IR data points. Finally, we interpolate the value at TI_m, compare it to both |M_m| and − |M_m| and select the polarity of M_m as the closest of the two points to the interpolated value at TI_m.

In a voxel containing tissue with heterogeneous T1-MD properties, the net signal can be described as a superposition of contributions from microscopic water pools with arbitrary T1 and MD properties. Assuming slow exchange between these subvoxel components, we can write the net signal attenuation as a function of the contrast weighting parameters b, TI, and TR:

where, R₁ = 1/T₁ and D̄ are the longitudinal magnetic relaxation rate and mean diffusivity, respectively, and p(R1,D̄) is the two-dimensional joint probability density function of R₁ and D̄ values within the voxel. A piecewise continuous approximation of p(R1,D̄) can be estimated numerically by solving the linear system M · p = y, with y and p representing column vectors containing the measured IR-IDE signal attenuations and the unknown signal fractions of T1-MD spectral components, respectively. The elements of the encoding matrix M, M_ij, are computed as the mean signal attenuations in the i-th experiment, S(b_i, TI_i, TR_i), evaluated over the 2D box function determined by the intervals [R_1,j, R_1,j+1] and [D̄j,D̄j+1] corresponding to the j-th spectral bin in the piecewise continuous approximation of p(R1,D̄).

Because this problem is poorly conditioned, additional regularization is required to obtain a stable solution. A straightforward approach is to solve the problem by least-squares minimization using L2-norm regularization and positivity constraints:

where, I is the identity matrix and λ is the L2-norm regularization parameter (Hansen, 1992).

RF field inhomogeneities can affect the efficiency of the adiabatic inversion yielding an inversion angle θinv≠180° and reducing the factor 2 in Equation (1) to (1 − cosθ_inv). In addition, dynamic processes such as magnetization transfer (MT) and cross-relaxation between immobile (bound) water associated with macromolecules and mobile (free) tissue water can produce a similar effect, further reducing the factor 2 in Equation (1) (Roscher et al., 1996; Gochberg et al., 1997; Does et al., 1998). Following the convention from Labadie et al. (2014), we can combine these processes together using the apparent inversion efficiency η defined as a percentage:

We propose a two-step process to fitting the IR-IDE data using Equation (4). First, we solve the constrained linear problem assuming perfect inversion (Equation 3) to derive an approximate solution. This step can be implemented using non-negative least-squares optimization, and the regularization parameter can be optimized with various methods (Canales-Rodríguez et al., 2021). However, since in this step we seek merely an approximate solution, it is not necessary to optimize the value of λ for each voxel. The approximate solution is then used as a starting point in solving the nonlinear optimization problem:

where, the elements of the encoding matrix M_η, are calculated as the mean signal attenuations in the i-th experiment, S_η(b_i, TI_i, TR_i), evaluated over the 2D box function determined by the intervals [R_1,j, R_1,j+1] and [D̄j,D̄j+1] corresponding to the j-th spectral bin in the piecewise continuous approximation of p(R1,D̄). This second step of the spectral reconstruction is implemented using constrained non-linear optimization. The final solution, [p^;η^] yields estimates for both the T1-MD spectral components and the parameter η.

From multiple measurements with different combinations of joint IR and IDE weightings, y, we can reconstruct the T1-MD spectrum, p^, in each voxel. These piecewise continuous spectra can approximate arbitrary probability density functions of T1-MD values that may occur in healthy and pathological tissues.

We performed Monte Carlo (MC) experiments to validate the numerical stability of the T1-MD spectral reconstruction and to empirically optimize an experimental protocol suitable to run on a clinical scanner. Considering the range of T1- and diffusion-encoding parameters achievable on a clinical scanner and the time constraints for scanning human subjects, we generated a protocol for acquiring 304 IR-IDE DWIs comprising all combinations of 19 T1-weightings, i.e., (TI,TR) pairs, and 16 diffusion weightings, i.e., b-values. We conducted MC experiments assuming η = 90% and using mixtures of normal or lognormal T1-MD distributions (Remin et al., 1989) as our ground-truth spectra. For each distribution, we computed the 304 IR-IDE ground-truth signal attenuations and we generated multiple instances of noisy data with signal-to-noise ratio (SNR) levels similar to those achievable in the human brain. From 500 independent instances of simulated noisy measurements, we reconstructed the corresponding T1-MD spectra, quantified the means and standard deviations in all spectral bins, and compared the results to the original ground-truth distributions. We also investigated the effect of the reconstruction grid size on the numerical stability of the spectral estimation and its ability to distinguish multiple peaks.

Next, we tested the imaging capabilities of our IR-IDE pulse sequence and assessed the accuracy of our T1-MD measurement (experimental design, contrast encoding, signal representation, and numerical estimation) in quantifying the physical properties of a well-calibrated homogeneous spherical diffusion MRI phantom that mimics the average relaxation properties of brain tissues. The phantom was constructed using a 40% polyvinylpyrrolidone (PVP) polymer solution (Pierpaoli et al., 2009; Wagner et al., 2017; Sarlls et al., 2018) with nominal diffusivity and T1 of 0.53 μm²/ms and 0.71 s, respectively, at room temperature. Diffusion within the phantom is expected to be uniform, isotropic and homogeneous, with a single mean diffusivity component and no signs of microscopic anisotropy or orientation dispersion. The homogeneity of the phantom at the macroscopic and microscopic scales allows us to concurrently assess potential sources of imaging artifacts such as ghosting, ringing, eddy current distortions, as well as potential errors in T1 and diffusion encodings due to source such as B1 inhomogeneities, gradient nonlinearities, concomitant fields, or eddy currents. We conducted an IR-IDE experiment in the phantom with the same protocol and reconstruction pipeline used for the in vivo exams. Finally, we compared the 2D T1-MD correlation spectroscopy images in the phantom to the nominal T1 and MD values.

We scanned three healthy volunteers using the IR-IDE protocol developed empirically based on the numerical simulation experiments. All subjects provided written and informed consent to participate in the study in accordance with a clinical protocol approved by the institutional review board (IRB) of the Intramural Research Program of the National Institute of Neurological Disorders and Stroke (NINDS). The MRIs were acquired on a conventional 3T Siemens Prisma scanner equipped with a 32-channel head RF coil and a gradient system that can achieve a maximum gradient strength of 80 mT/m/axis. We acquired 304 IR-IDE DWIs in 19 separate scans with different T1-weightings: one scan with no-IR and TR = 12 s, and 18 scans with TI values between 0.05 and 5 s and TR values between 1.63 and 10 s. For each scan, we first played out three discarded acquisitions (3 TRs) to establish the steady-state longitudinal magnetization followed by 16 consecutive acquisitions during the following TRs with a range of IDE weightings b = 50–3,600 s/mm². All volumes were acquired with single-shot multislice EPI readout with TE = 98 ms, a field-of-view (FOV) of 220 mm on an 88 × 88 imaging matrix, a GRAPPA factor of 2, no partial Fourier sampling, 2.5 mm in-plane resolution, 20 axial slices with a 5 mm slice thickness grouped in two concatenations. The non-cubic voxel shape 2.5 × 2.5 × 5.0 mm was chosen to mitigate imaging artifacts due to Gibbs ringing while maintaining a high voxel SNR. The total scan time of the IR-IDE protocol was 51 min. All volunteers were instructed to minimize motion during the entire duration of the experiment. Additional scans were also acquired: a 30-direction DTI scan with b = 1,100 s/mm², and the same EPI acquisition parameters as the IR-IDE scans; a 1-mm sagittal MP-RAGE (TE/TI/TR = 3.3/1,100/2,530 ms); and a 1-mm T2-weighted (T2W) Turbo-spin echo image (TE/TR = 72/8,000 ms) to serve as an anatomical reference for EPI distortion correction.

We processed all diffusion MRI data sets using the TORTOISE software package to correct for EPI distortions due to eddy currents and magnetic field inhomogeneities (Pierpaoli et al., 2010) and to register the DWIs to the corresponding T2W structural reference scan. We estimated tissue SNR levels by dividing the baseline signals (no IR, b = 0.05 ms/μm²) averaged in tissue-specific regions-of-interest by the signal standard deviation measured in a noise region outside the brain. From the corrected IR-DWIs we computed voxelwise T1-MD spectra using 144 spectral components defined by bins with 12 MD values logarithmically spaced between 0.3 and 3.0 μm²/ms, and 12 T1 values logarithmically spaced between 0.25 and 3.3 s. We repeated the analysis using different grid sizes with logarithmic spacing in order to better visualize signal components in specific spectral bands. We derived marginal distributions of subvoxel T1 values and MDs by integrating the T1-MD spectra along each dimension and compared the variability of the estimated spectra across subjects. Furthermore, based on visual inspection of whole-brain T1-MD distributions we manually selected spectral components that are consistent with different tissue types and integrated over the corresponding spectral bins to determine the signal fractions of these components in all subjects.

MC simulation results suggest that we can estimate subvoxel multicomponent T1-MD spectra from a sufficiently large number of IR-DWIs with SNR levels that can be achieved on clinical MRI scanners. Figures showing examples of MC simulations using our clinical protocol and ground-truth spectra with one, two, and three components similar to those observed in vivo at different SNR levels are included as Supplementary Material. Figure 2 shows MC simulation results for a ground-truth spectrum containing a mixture of three components. The apparent inversion efficiency can be measured with very good precision and accuracy. The peak locations and relative signal fractions in the mean estimated and ground-truth correlation spectra and in the corresponding 1D marginal distributions are generally in good agreement. At higher SNR levels, the uncertainties in estimating the peak locations and amplitudes are lower, reducing the variance of the reconstructed spectra and the blurring in the estimated mean spectrum.

Apart from noise, the accuracy of the spectral reconstruction in our MC simulations is also affected by the finite number of measurements and the use of constraints in solving the poorly-conditioned problem. The constraints in Equation (5) make it difficult to establish analytical relations between the spectral resolution, accuracy, or precision and factors such as measurement noise or encoding parameters (b, TI, and TR) that apply generally to any ground-truth distribution. Nevertheless, the performance of the spectral reconstruction is relatively unaffected by our choice of the reconstruction grid (Figure 3). Because the estimated normalized spectra are piecewise continuous probability density functions (i.e., 2D box functions), the reconstruction grid (i.e., spectral bins) does not affect the characteristics of the solution. While in this study we used an empirical approach based on the MC experiments to derive a protocol suitable for use on a clinical scanner, the experimental design may be improved in the future with advanced sampling methods and optimization algorithms.

Images obtained in the PVP phantom showed single-peak T1-MD spectra throughout the phantom (Figure 4), with MD ≈ 0.55μm²/ms and T1 ≈ 0.70s, in good agreement with the nominal values (Pierpaoli et al., 2009; Wagner et al., 2017; Sarlls et al., 2018). These results illustrate the ability of our imaging method to accurately encode and map T1-MD properties across large fields-of-view. The accuracy in estimating the MD of the PVP phantom confirms that the diffusion gradient amplitudes are scaled accurately to obtain the desired b-values. Meanwhile, the relatively small apparent inversion efficiency values suggest that RF inversion pulses are accurately calibrated and the T1-encoding in our sequence is not significantly affected by B1 inhomogeneities. Slight variations in the spatial-spectral uniformity of the estimated spectra may be caused by imaging distortions and unaccounted diffusion encoding contributions from magnetic field inhomogeneities, eddy currents, concomitant fields, or gradient non-linearities.

IR-IDE DWIs showed good SNR in all study participants. Typical SNR values computed from non-attenuated signals were 250 in gray matter (GM), 125 in WM, 105 in the basal ganglia (BG), and 370 in cerebrospinal fluid (CSF). The relatively high SNR in GM is likely due to partial volume contributions from CSF, while the use of parallel imaging may also bias these tissue SNR estimates. Throughout the brain, signal attenuations due to diffusion and T1 contrasts showed a large dynamic range which was further increased after phase-sensitive IR correction (Figure 5).

We found predominantly T1-MD spectra with one, two, or three distinct components in the brains of healthy volunteers. Figure 6 shows the raw data points and reconstructed normalized T1-MD correlation spectrum in a representative brain tissue voxel along with the residuals to the spectral reconstruction as a function of the encoding parameters. The independence of the residuals on the T1 and diffusion encoding suggests that the spectral representation explains the data well. In vivo estimates of the inversion efficiency parameter, η, showed values of ~ 89% in WM, 95% in GM, and 98% in CSF (Figure 7), in good agreement with an earlier study (Labadie et al., 2014). Maps of normalized T1-MD correlation spectra and corresponding 1D marginal distributions revealed signal components that were spatially consistent with specific tissue types such as WM, GM, BG, and CSF (Figure 7).

The CSF signal showed a single peak (Figure 7, blue contour) with long T1 and large MD values (Kwong et al., 1991) that was well-isolated from parenchymal spectral components. Due to the relatively coarse spatial resolution, many voxels in the cortex and around the ventricles showed both CSF and parenchymal components. These partial volume contributions could be promptly separated in the reconstructed spectra. Even though the CSF peak can be easily distinguished in both spectral dimensions, it has a narrower range of MD values centered around ≈ 3.0μm²/ms. The larger range of the CSF T1 values 2.0–3.3 s may reflect partial volume variations across measurements with different TIs due to subject/physiological motions or the inflow of fresh blood (i.e., non-inverted spins) from neighboring slices especially in GM (Trampel et al., 2004) during the T1 encoding duration. Diffusion encoding on the other hand occurs during a fixed short duration of the evolution of transverse magnetization and is localized to the excited slice.

In the parenchyma, WM showed the largest heterogeneity of subvoxel T1-MD properties, revealing two distinct peaks. The major peak (Figure 7, green outline) was relatively broad with intermediate values of T1 1.0–1.3 s and MD 0.85–1.00 μm²/ms. The second peak with shorter T1 values ≈ 250 ms was specific to the major WM pathways and accounted for ≈ 10.2% of the WM signal (Figure 7, magenta outline). It may reflect processes associated with myelin water as suggested in previous in vivo studies (Labadie et al., 2014). The MD properties of the short-T1 peak were similar to those of other parenchymal components suggesting that this peak may be an indirect measure of myelin water via MT or chemical exchange (Avram et al., 2010). The short-T1 WM component could be better visualized using a higher resolution reconstruction grid in Figure 8. The reconstructed spectra in WM, and more generally in the brain parenchyma, did not reveal very low MD components that may suggest the presence of a very restricted diffusion compartment (i.e., dot compartment; Tax et al., 2020).

In GM we observed a single component, with significantly longer T1 values 1.30–2.08 s and slightly lower MD values compared to WM (Figure 7, red contour). The lowest parenchymal MD values 0.56–0.69 μm²/ms were specific to the BG, including the putamen and caudate nucleus (Figure 7, yellow contour), in agreement with previous in vivo studies (Avram et al., 2019). The localization of low MD components in the BG is evident in the marginal MD distribution components (Figure 7, left column). The BG also contained spectral components found in other tissues.

The one-dimensional (1D) marginal distributions derived from the T1-MD correlation spectra were in very good agreement with results from other T1 (Labadie et al., 2014) and MD (Avram et al., 2019) relaxation spectroscopic MRI studies, further supporting the accuracy of the T1 and MD preparations and the spectral reconstruction. The signal fractions obtained by integrating the manually selected spectral regions-of-interest (ROIs) corresponding to 1 CSF and 4 parenchymal components were consistent across the three study participants (Figure 9) and showed slight improvements in anatomical specificity compared to the 1D marginal distributions (Figures 7, 8). The ad hoc approach for defining the spectral ROI components in Figure 7 will be improved and automated in future studies with larger cohorts of subjects. While these preliminary findings did not reveal any remarkable T1-MD correlations in the healthy brain, they encourage further studies with larger cohorts of healthy volunteers and patients to assess the reproducibility and establish the clinical applicability and utility of T1-MD CS-MRI.

The spatial-spectral characteristics of T1-MD distributions did not depend on the locations or sizes of the spectral bins in the reconstruction (Figures 7, 8), as expected from our numerical simulation results (Figure 3). Estimating T1-MD spectra with different reconstruction grids can improve visualization of specific spectral bands. For example, Figure 8 reveals that the T1 values of the short-T1 component in WM are in the range of 0.13–0.35 s. In general, reconstructions using denser grids require longer processing times and can be more sensitive to measurement noise as they aim to estimate more unknowns (i.e., spectral components). In addition, reconstructions using dense grids distribute the spatial-spectral information over a large number of spectral components maps, making it more difficult to identify spatially coherent anatomical features. Future studies will use numerical simulations to optimize not only the reconstruction grids but also the experimental protocols (TI-b weightings) for specific clinical applications based on high-quality T1-MD CS-MRI data sets measured in larger cohorts of healthy volunteers and patients.