Fast and reproducible in vivo T1 mapping of the human cervical spinal cord

To develop a fast and robust method for measuring T1 in the whole cervical spinal cord in vivo, and to assess its reproducibility.


INTRODUCTION
The longitudinal relaxation time T 1 is one of the most fundamental quantitative parameters in MRI. Several studies have investigated the biological correlates of T 1 . It is well-established that T 1 is dependent on myelin content (1). Relaxation of myelin water occurs faster than nonmyelin water; therefore, tissues with higher myelin content show a lower average T 1 . However, T 1 has been shown to depend also on additional microstructural features with a constant total myelin volume, such as water content (2), axonal size in white matter (3), and iron concentration in gray matter (4).
Despite the lack of specificity to a single particular biological feature, T 1 is sensitive to changes in tissue microstructure caused by pathologies and inflammatory events. It also provides a quantitative measure comparable across subjects and centers, which is more informative than visual examination of conventional T 1weighted images. Precise and robust characterization of T 1 in vivo is therefore of great importance. Moreover, the accurate knowledge of T 1 serves as the basis for other quantitative MR methods (e.g., perfusion and quantitative magnetization transfer imaging) and in the optimization of imaging sequence parameters.
Several techniques are available to estimate T 1 in vivo, especially in the brain. Considerable progress has been made in the development of fast mapping techniques: methods such as the Look-Locker inversion recovery (IR) (5,6) and the variable flip angle (VFA) (7,8) are particularly appealing for their ability to achieve T 1 estimates with large coverage in a few minutes of scan time.
However, great variability among T 1 estimates can be found in literature, which has been ascribed to the sensitivity of the various methods used to site-specific factors such as radiofrequency (RF) field uniformity, RF pulse imperfections, and incomplete spoiling (9). Thus, the time-consuming IR T 1 mapping approach is still regarded as the reference ("gold standard") method for in vivo T 1 measurements.
Development of quantitative MRI techniques for the spinal cord, including the characterization of T 1 , generally lags behind the brain (10,11), despite the crucial involvement of the spinal cord in several diseases, such as multiple sclerosis (12)(13)(14), amyotrophic lateral sclerosis (15), spinal cord injury (16), and neuromyelitis optica (17). Current techniques for measuring T 1 in the spinal cord rely on the mere application of brain protocols to the cervical level of the spinal cord, mostly using the VFA method (18). However, the VFA method makes use of 3D (or 2D) spoiled gradient-echo acquisitions, which are inherently more sensitive to intrascan motion artifacts when compared with other approaches such as singleshot echo-planar imaging (EPI). If not properly addressed, these artifacts could propagate into the T 1 maps, especially in spinal cord applications, in which effects of physiological noise and subject motion are exacerbated. Additionally, several volumes (with varied flip angle) should be acquired to avoid noise bias in the fitting (19), complete magnetization spoiling has to be ensured (20), and an accurate map of the transmitted RF field (B 1 ) has to be obtained to correct for the actual flip angle at the voxel level (21). These facts contribute to preventing the design of fast, accurate, and robust protocols for T 1 mapping in the spinal cord using the VFA method.
In contrast, IR-based methods are less sensitive to the aforementioned limitations (9), and thus represent a more robust approach to for T 1 mapping in challenging environments such as the spinal cord, where multiple noise sources and field inhomogeneities coexist.
In this study, we introduce an IR-based T 1 mapping method for the spinal cord, making use of the sliceshuffling scheme initially developed for application to the brain (22,23). Our method allows for in vivo T 1 mapping of the whole cervical spinal cord in a clinically acceptable scan time. We demonstrate its intra-and intersubject reproducibility in a small cohort of healthy subjects.

METHODS
The slice-shuffling mechanism consists of cycling the order in which different slices are acquired over sequence repetitions. Here, we used the slice-shuffling mechanism in the context of T 1 mapping (22,23).
An IR sequence was developed for T 1 measurements in the spinal cord, making use of zonally oblique-magnified multislice (ZOOM) EPI (24)(25)(26), which allows artifact-free multislice imaging of small structures using a single-shot EPI readout. Slices are acquired with an interleaved order, allowing a time interval between contiguous slices excitation (TR) long enough for restoration of longitudinal magnetization after each oblique spin-echo pulse pair. This results in N p groups (i.e., packages) of N spp ¼ N s /N p maximally spaced-out slices acquired every TR, in which N s is the total number of prescribed slices (Fig. 1a). An IR experiment is implemented using nonselective adiabatic inversion pulses, applied before the acquisition of each FIG. 1. a: Schematic of ZOOM-EPI multislice acquisition in the whole cervical cord, from levels C1 to C7. Slices are grouped in packages (identified by colors); a package is acquired per each TR. Inversion-recovery ZOOM-EPI is implemented by adding a global inversion pulse before each package acquisition (i.e., in each TR). b: Between acquisitions of different packages, a recovery time (T rec > 5 s) is added to ensure full recovery of magnetization. Slice order is shuffled within packages to avoid TR lengthening in a multi-TI experiment, producing a number of effective TIs equal to the number of slices per package (N spp ).
package. Hyperbolic secant adiabatic inversion is used to achieve robust inversion over the target volume against B 1 and B 0 inhomogeneities (27).
Multiple inversion time (TI) data are acquired by shuffling the slice-acquisition order within the package for any of the M given delays between the inversion pulse and the first excited slice (TI app ) (i.e., for a given TI app , a total of N spp effective TIs are obtained). To avoid contamination between tilted slice excitations and inversion pulses, the TR is extended by a recovery time T rec ¼ 6 s (Fig. 1b).
The combination of IR preparation with a rapid ZOOM-EPI readout allows the exploitation of the sliceshuffling mechanism to reduce dead time when large slice coverage is required. For a given stack of N s slices grouped in N p packages, and M different TI app , resulting in N spp Â M different effective TIs (TI eff ), the total protocol time T prot is where uniform temporal slice spacing Dt s is assumed. The sampled TIs can be chosen by rearranging Dt s , TI app , and potentially M, according to the specific applications that dictate coverage N s ¼ N p Â N spp , and scan time available T prot . Although a theoretical investigation of protocol optimization is beyond the scope of this article, Equation [1] shows how different settings of N spp , N p , and M for the same N s can be used to reduce T prot , therefore providing the basis to exploit the interplay among sequence parameters N s , N p , N spp , and sampling scheme parameters M, TI app , and Dt s in optimization procedures. Five healthy volunteers (age range 27-37, 4 males, 1 female) were enrolled in the study. Imaging was performed on a Philips Achieva 3T (Philips Healthcare, Best, the Netherlands) MRI system with a 16-channel neurovascular coil, using parallel transmission technology.
The whole cervical spinal cord (i.e., C1-C7) was imaged using the following MR parameters: field of view (FOV) ¼ 64 Â 48 mm 2 , in-plane voxel size ¼ 1 Â 1 mm 2 , EPI train length ¼ 63, echo time (TE) ¼ 22 ms, partial Fourier factor ¼ 0.6, 24 slices, 5 mm thick. Magnetization recovery was sampled using the IR-ZOOM-EPI sequence described previously, at 12 TIs obtained by grouping slices in N p ¼ 4 packages and using M ¼ 2 different TI app of 100 and 1300 ms, and slice temporal spacing Dt s ¼ 200 ms, to produce the following TI eff ¼ 100, 300, 500, 700, 900, 1100, 1300, 1500, 1700, 1900, 2100, and 2300 ms. A "noise-only" scan obtained without RF and gradients was added to characterize noise standard deviation (r) in the quantification. Total protocol duration was 7 min, 6 s. Validation of the protocol against a standard single slice IR in phantoms is shown in Supporting Figures S1 and S2 and Supporting Table S1.
The protocol was repeated in a separate session for each subject to assess the reproducibility. Within-subject motion across time points was corrected slice-wise by means of 2D linear transformations with 3 degrees of freedom (28,29) using a model-based image registration approach (30).
Straightening of spinal cords was performed on IR-ZOOM-EPI and 3D-GRE data based on (31), enabling inherent coregistration among different modalities and facilitating registration of a spinal cord template to 3D-GRE volumes in subsequent analysis.
The following mono-exponential signal-recovery model was fitted to magnitude data using maximum likelihood estimation assuming Rician distributed noise, to estimate the equilibrium magnetization M 0 and T 1 . s was obtained from the noise-only scan after smoothing image intensities with a moving average filter (3 Â 3 kernel size) and correcting for noise floor bias (32) as follows: where g is the image intensity of the noise-only scan, after smoothing. Spinal cord ROIs were defined automatically via registration of a spinal cord template (23,24) to each individual 3D-GRE scan, separately for each session. Template registration was performed using the spinal cord toolbox (25), and refined by supplying cord white matter (WM) and gray matter (GM) masks obtained using the method found in (33). Four different ROIs were selected: (i) GM obtained from segmentation, (ii) WM dorsal column obtained by merging the fasciculus cuneatus and fasciculus gracilis atlases from the template, (iii) WM left lateral column (LC L ), and (iv) WM right lateral column (LC R ) obtained by merging the respective spinothalamic, spinoreticular, rubrospinal, and lateral corticospinal tracts atlases from the template.
The mean T 1 and standard deviation were measured in each ROI, and in the whole cord, for both sessions. The intersession coefficient of variation (COV) was calculated for each ROI using where T i;j 1 refers to the mean T 1 of subject i (with i ¼ 1,2,..,5) for session j (j ¼ 1,2) in a specific ROI.
Reproducibility of T 1 estimates at the regional level was assessed using Bland-Altman plots, linear regression, and correlation analysis.
The intraclass correlation coefficient (ICC1,2) was calculated for ROIs pooled among subjects using a two-way random effects model (34,35). The ICC measures the fraction of the total variability attributable to differences among subject-specific ROIs.

RESULTS
An example of the acquired IR data is shown in Figure 2 for different cervical levels.
T 1 maps at different vertebral levels of the cervical spinal cord (from C1 to C7) are shown in Figure 3 in a single subject. T 1 values appear very homogenous across all of the cord levels and show very little variability between scan and rescan maps.
Whole-cord T 1 mean, averaged across sessions and subjects, was 1142.3 (6148.2) ms. Regional T 1 values across subjects and sessions are given in Table 1. The average T 1 was 1108.5 (677.2) ms for LC L , 1110.1 (683.2) ms for LC R , 1150.4 (6102.6) ms for dorsal column, and 1136.4 (690.8) ms for GM. Intersession COV for different ROIs was 0.93% for LC L , 0.89% for LC R , 0.74% for dorsal column, and 1.02% for GM. The COV for the whole-cord ROI was 0.95%.
Regional T 1 estimates showed good scan-rescan correlation, as demonstrated in Figure 3a. The Pearson correlation coefficient was 0.89 (P value < 0.01), and the slope of linear regression was 0.73 (confidence interval 0.54-0.92). The Bland-Altman plot in Figure 3b shows a negligible bias among T 1 estimates (2 ms), and very narrow 95% limits of agreement (approximatively 6 20 ms). The intraclass correlation coefficient for regional T 1 was 0.88 (P value < 0.01, confidence interval 0.71-0.95).

DISCUSSION AND CONCLUSIONS
In this study, we demonstrated a new, simple, fast, and reproducible method for mapping T 1 in the whole cervical spinal cord in vivo without deviating from the standard IR approach. Quantitative characterization of spinal cord microstructure is important in a variety of neurological disorders. More effort is required in developing robust methods to assess it, as many aspects of spinal cord microstructure are as yet uncharacterized in vivo with the current state of the art of quantitative MRI.
The proposed sequence, IR-ZOOM-EPI, combines advantageous solutions for T 1 mapping in the spinal cord. The use of a ZOOM-EPI readout to perform reduced FOV single-shot acquisition, which has previously been applied at 3 T in diffusion studies (36,37), while limiting distortions and freezing intrascan motion, is essential to achieve clinically feasible protocol durations. In fact, a previous in vivo T 1 mapping study in the spinal cord using IR required 20 min per slice and used a large FOV (18). In addition, the ZOOM-EPI readout is compatible with the magnetization preparation of an ideal IR experiment, in which inversion is performed by a nonspatially selective pulse. Because a TR constraint (TR > > T 1 ) is inherently in place with ZOOM-EPI to avoid cross-contamination from the oblique excitation pulse, the use of such an inversion pulse has a limited time penalty compared with the normal spin-echo sequence (the maximum TR increment is always lower than the maximum TI eff chosen). The adiabatic inversion is also robust against magnetic field inhomogeneities, which are exacerbated in the spinal cord (11), and effects of the imperfect slice profile of the inversion pulse, thought to contribute to the variability of T 1 measurements (9).
Together with the adiabatic inversion, slice shuffling enables the optimal use of scan time (22,23,38), as a  range of effective TIs can be produced for the same nominal TR.
The combination of these factors into a single sequence allows the design of time-efficient T 1 -mapping IR protocols for the spinal cord, of comparable duration to those using the VFA method with linear approximation fitting.
Quantitative T 1 values measured here lie within the expected range for tissues in the central nervous system at 3 T (39,40), with WM values being on the upper side of the reported values. This could be explained by the presence of larger WM axons (41), which leads to higher T 1 values for constant myelin volume fraction (3), or the presence of myelinated gray matter, as shown by histological findings (12,42). Other spinal cord studies, in which the T 1 was not the primary parameter of investigation, have indeed reported even higher values than those obtained here (43)(44)(45).
From the scan-rescan experiments we obtained excellent measures of reproducibility. Intersession COV was approximately 1% for all of the ROIs considered. This suggests that the T 1 can be measured precisely in the spinal cord using a relatively short protocol. Moreover, the high ICC values suggest that most of the variabilities in the T 1 estimates among subjects and ROIs derive from biological differences rather than measurement errors.

Limitations
Contrary to brain T 1 studies, we did not observe consistent contrast between WM and GM. As previously reported (18), this may be the result of the axial resolution used in this study (1 Â 1 mm 2 ), which may not allow for clear delineation of the GM/WM boundaries, as well as contamination from cerebrospinal fluid flow. The use of higher resolution is warranted for pursuing detailed tissue specific analysis. Use of a 5-mm slice thickness is common in axial spinal cord imaging, and justified by the smooth variation of spinal cord anatomy in the cranial-caudal direction. However, the use of thinner slices could be beneficial for improving the accuracy of the tissue-specific characterization, as well as for the efficacy of motion correction techniques. Residual misalignments can in fact further confound the resultant WM/GM separation in parametric maps.
The periodic alternation of vertebrae and intervertebral disks along the spinal cord produces intensity signal variations in the slice direction (as seen along the columns of Fig. 2). The trend is evident also on maps displayed in Figure 3, with T 1 slightly decreasing with the cervical level, which is consistent with recent findings at ultrahigh field (46). The implementation of techniques to mitigate signal intensity variation along the rostro-caudal direction, such as slice-wise z-shimming (47), should be considered in future studies to precisely assess the impact of this effect on T 1 estimation.

Future Improvements and Conclusions
The approach proposed here can be used in multimodal studies to better characterize spinal cord microstructure. IR-ZOOM-EPI has the potential to be extended to other spinal cord levels, as ZOOM-EPI has proven successful in diffusion studies at the lumbar level (48).
Additionally, more room for optimization is available in terms of parameter precision, via accurate selection of TIs with protocol optimization techniques (49,50) and acquisition time through combination with ultrafast imaging techniques (e.g., simultaneous multislice imaging (51)). The singular interplay among sequence parameters already allows scan time to be gained without needing to use acceleration techniques (i.e., N spp and M can be arranged to provide a time-efficient protocol even when higher N s are used). Formal optimization that accounts for protocol duration and TI selection could be devised starting with Equation [1], to produce an optimal T 1 protocol in the spinal cord, simultaneously addressing the issues of spatial coverage, SNR, protocol time, and T 1 estimate precision. Additional time gains could be used, in turn, to increase SNR through signal averaging for higher resolution imaging.
Alternative emerging approaches for T 1 mapping, such as MR fingerprinting (52), could allow further scan-time reduction, although its implementation to the spinal cord has not yet been proposed, and only applications to the brain are currently available.
In conclusion, reproducible measurements of T 1 can be obtained in the whole cervical spinal cord with a rapid protocol. Our method enables more comprehensive T 1 mapping studies in the spinal cord for different pathologies.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article. Fig. S1. Comparison of T 1 estimates in phantoms between standard IR (a) and IR-ZOOM-EPI (b). The IR-ZOOM-EPI is repeated separately on each phantom (with the respective FOV highlighted in the corner of each map). N s 5 6 slices were acquired. The central slice is shown here for the purpose of comparison. Fig. S2. Correlation plots between T 1 estimates from standard IR and IR-ZOOM-EPI with T rec 5 10 s (a), T rec 5 8 s (b), and T rec 5 6 s (c). In each plot, the widths of the marker in each direction (horizontal and vertical) correspond to the width of the relative distribution (standard IR and IR-ZOOM-EPI) within the 1st and 99th percentile. The identity line (indicating an ideal linear relationship), together with parameters for the linear fit of T 1 from IR-ZOOM-EPI against T 1 from IR single slice, are shown. The bar graph (d) shows a quantitative comparison between T 1 estimates in all of the phantoms. The bar height represents the mean T 1 value, whereas the error bars indicate the 99th percentile of the distribution. No differences between standard IR and all versions of IR-ZOOM-EPI were detected by a paired ttest. Table S1. Mean and Standard Deviation of T 1 Estimates in the Phantoms for Standard IR Single Slice, IR-ZOOM-EPI With T rec 5 10 s, IR-ZOOM-EPI With T rec 5 8 s, and IR-ZOOM-EPI With T rec 5 6 s