Multi-parametric quantitative spinal cord MRI with unified signal readout and image denoising

Multi-parametric quantitative MRI (qMRI) of the spinal cord is a promising non-invasive tool to probe early microstructural damage in neurological disorders. It is usually performed by combining acquisitions with multiple signal readouts, which exhibit different thermal noise levels, geometrical distortions and susceptibility to physiological noise. This ultimately hinders joint multi-contrast modelling and makes the geometric correspondence of parametric maps challenging. We propose an approach to overcome these limitations, by implementing state-of-the-art microstructural MRI of the spinal cord with a unified signal readout. We base our acquisition on single-shot echo planar imaging with reduced field-of-view, and obtain data from two different vendors (vendor 1: Philips Achieva; vendor 2: Siemens Prisma). Importantly, the unified acquisition allows us to compare signal and noise across contrasts, thus enabling overall quality enhancement via Marchenko-Pastur (MP) Principal Component Analysis (PCA) denoising. MP-PCA is a recent method relying on redundant acquisitions, i.e. such that the number of measurements is much larger than the number of informative principal components. Here we used in vivo and synthetic data to test whether a unified readout enables more efficient denoising of less redundant acquisitions, since these can be denoised jointly with more redundant ones. We demonstrate that a unified readout provides robust multi-parametric maps, including diffusion and kurtosis tensors from diffusion MRI, myelin metrics from two-pool magnetisation transfer, and T1 and T2 from relaxometry. Moreover, we show that MP-PCA improves the quality of our multi-contrast acquisitions, since it reduces the coefficient of variation (i.e. variability) by up to 15% for mean kurtosis, 8% for bound pool fraction (BPF, myelin-sensitive), and 13% for T1, while enabling more efficient denoising of modalities limited in redundancy (e.g. relaxometry). In conclusion, multi-parametric spinal cord qMRI with unified readout is feasible and provides robust microstructural metrics with matched resolution and distortions, whose quality benefits from MP-PCA denoising, a useful pre-processing tool for spinal cord MRI.


Introduction
The spinal cord is a small but functionally important structure of the human central nervous system, affected in several common disorders. These are often associated with high disability (Hendrix et al., 2015), and include: multiple sclerosis (Ciccarelli et al., 2019), amyotrophic lateral sclerosis (van Es et al., 2017), spinal cord injury (Ahuja et al., 2017) and many others (Lorenzi et al., 2019). Routine anatomical magnetic resonance imaging (MRI) plays an important role in the diagnosis and management of these conditions (Kearney et al., 2015). However, it only offers macroscopic descriptors of tissue damage that lack specificity for pathophysiology, have limited prognostic value and fail to guide treatment and rehabilitation personalisation (Cohen-Adad, 2018;Stroman et al., 2014;Wheeler-Kingshott et al., 2014). The gradual adoption of quantitative MRI (qMRI) techniques may help overcome the limitations of conventional anatomical MRI. Based on either well-validated biophysical models or parsimonious signal representations , qMRI promises to provide estimates of biologically meaningful characteristics, which would make parametric maps vendor-independent (Cercignani and Bouyagoub, 2018). The latest multimodal qMRI techniques exploit the complementary information from different contrasts (De Santis et al., 2016;Stikov et al., 2015), for example relaxometry and diffusion, to better quantify the parameters of tissue microstructure (Lemberskiy et al., 2018;Ning et al., 2019;Slator et al., 2019;Veraart et al., 2018). qMRI of the spinal cord is increasingly popular (Battiston et al., 2018a;Battiston et al., 2018b;By et al., 2017By et al., , 2018Duval et al., 2017;Grussu et al., 2019;Grussu et al., 2015;Ljungberg et al., 2017;Massire et al., 2016;Schilling et al., 2019;Taso et al., 2016) due to recent advancements in scanner hardware (Barry et al., 2018;Duval et al., 2015) and analysis software (De Leener et al., 2017). However, its development is currently hampered by the following two challenges.
Firstly, multi-contrast qMRI in the spinal cord typically relies on specialised techniques with dedicated signal readout for each contrast (Duval et al., 2017;Massire et al., 2016;Taso et al., 2016). The variety of readouts is not compatible with joint computational modelling of voxel-wise multi-contrast signals, and also limits the alignment of multimodal metrics due to different distortions and susceptibility to physiological noise (Campbell et al., 2018). The second major challenge is related to the fact that data quality in spinal cord imaging remains lower compared to the brain. This is due to the need for high spatial resolution (the spinal cord cross sectional area is about 1 cm 2 ), which is challenged by artifacts from pulsation (Morozov et al., 2018;Summers et al., 2006) and local magnetic field inhomogeneities Verma and Cohen-Adad, 2014). Improving the intrinsic quality of qMRI data is therefore imperative to facilitate the application of the latest qMRI techniques to the spinal cord, which are still in their infancy as compared to those in the brain (Cohen-Adad, 2018;Wheeler-Kingshott et al., 2014).
In this paper, we propose a unified acquisition for state-of-the-art multimodal qMRI of the spinal cord that addresses both challenges. Our protocol relies on a unified signal readout based on single-shot spin echo planar imaging (EPI) with reduced field-of-view (rFOV), which provides multi-contrast images, each having the same intrinsic resolution and susceptibility artifacts (i.e. distortions). Importantly, the unified acquisition enforces the same noise statistics across multiple signal contrasts, thus enabling overall data quality enhancement via Marchenko-Pastur (MP) Principal Component Analysis (PCA) denoising (Veraart et al., 2016a;Veraart et al., 2016b). MP-PCA is designed to denoise redundant qMRI data, for which the number of measurements within a given patch of voxels is much larger than the number of linearly-independent contributions to the signal from this patch (i.e., the number of principal components in the low-rank signal representation). MP-PCA applies naturally to acquisition-rich qMRI techniques, such as diffusion MRI, but its application in relaxometry is more challenging due to the limited number of measurements. Here, we exploit synergies between qMRI contrasts enabled by the unified readout, and test whether the combination of less redundant acquisitions with more redundant ones improve MP-PCA denoising of the former.
Our results demonstrate that a unified acquisition enables reliable multicontrast characterisation of spinal cord microstructure. A unified readout provides metrics that offer complementary information and are intrinsically co-aligned and matched for distortions, without significant losses of resolution for the least redundant modalities. Importantly, data from multiple vendors also demonstrate that MP-PCA is an important pre-processing tool for spinal cord qMRI, with the potential of bringing qMRI one step closer to the clinic by improving the repeatability of microstructural metrics. A unified readout practically increases the redundancy of multi-contrast qMRI data sets, enabling efficient denoising across a range of MRI contrasts.

Background on MP-PCA
This section outlines key concepts behind MP-PCA denoising required for the understanding of the remainder of the paper.
MP-PCA denoises noisy input matrices = [ , ] of size × constructed by arranging measurements along rows from neighbouring voxels along columns, such that < without loss of generality. MP-PCA requires redundant qMRI protocols, i.e. such that ≫ , with being the number of underlying linearlyindependent signal sources. In the absence of noise, identifies the number of nonzero singular values of , which is much smaller than when is redundant. This is the standard assumption behind the PCA-based dimensionality reduction, where one is looking for a P-dimensional hyperplane to represent an M-dimensional measurement. Importantly, in the presence of noise the P-dimensional hyperplane extends in all remaining dimensions, such that becomes full-rank (Veraart et al., 2016a). The denoising problem is then equivalent to identifying the original hyperplane and its dimensionality .
MP-PCA employs random matrix theory to identify information-carrying components within the { | = 1, … , } singular values of , and to draw a threshold between pure-noise and signal-carrying singular values (Johnstone, 2006). As the pure-noise principal components in the limit ≫ 1 are distributed according to the universal MP distribution, the denoising algorithm identifies the MP distribution within the eigen-spectrum of (Veraart et al., 2016b), and then sets the corresponding − singular values to zero. The information-carrying above the right edge of the MP distribution constitute the denoised signal. The output of MP-PCA is den = [ , den ], a denoised version of , as well as the simultaneously estimated noise standard deviation and the number of signal components (or signal generators) . Note that the estimated will be in general smaller than the original hyperplane dimensionality (without the noise), since some of the informative principal components may fall into the MP noise bulk (Johnstone, 2006). However, we point out that there are no noisefree measurements, and the estimated represents the effective dimensionality of the informative part of the signal, which remains above the noise level. For details of application of random matrix theory and MP-PCA denoising in MRI, please see (Veraart et al., 2016a;Veraart et al., 2016b).
To validate this (or any) denoising method, one can study the distribution ( , ) of the normalised residuals , = ( , den − , ) . (1) A perfect denoising method removes only Gaussian noise and does not remove tissue anatomy. Hence, as a way to monitor denoising quality, we will check that MP-PCA residuals are Gaussian, i.e.
such that the histogram ( , ) plotted against , 2 in the semi-log scale is a straight line with slope − 1 2 (Veraart et al., 2016b). If tissue anatomy is removed, the histogram plotted as such will not be a straight line, typically blowing up at the tail. Importantly, the line will have a steeper slope if not all noise is removed assuming the no tissue anatomy is spoiled and that is correctly estimated, i.e. (the removed noise distribution is narrower, and has a smaller variance than the one actually present in the measurement). The performance of denoising increase as the slope approaches − 1 2 from below (Veraart et al., 2016b), implying that a larger fraction of noise-induced signal variability is mitigated.

Methods
We synthesised multimodal MRI scans encompassing modalities with different redundancy, emphasising protocols that could be realistically implemented in the spinal cord in vivo, and evaluated the performance of MP-PCA denoising when performed on each modality independent or on multiple modalities jointly.
We also acquired multi-contrast MRI data with unified readout on scanners from two vendors (vendor 1: Philips; vendor 2: Siemens), and characterised the quality of several qMRI metrics as obtained following MP-PCA denoising or without denoising. The shared readout enables the assessment of whether denoising modalities characterised by limited redundancy can be improved if these are denoised jointly with more redundant acquisitions.
In the following sections, we will describe simulations first, as these provide the context for the interpretation of findings in vivo. All analyses were performed using inhouse scripts, which are made openly available (http://github.com/fragrussu/PaperScripts/tree/master/sc_unirea dout).
Firstly, we used NiftyReg (http://niftyreg.sf.net) reg_resample (Modat et al., 2010) with default options to downsample the voxel-wise volume fractions WM , GM and CSF to a resolution that is plausible for quantitative MRI of the spinal cord based on EPI Duval et al., 2015;Grussu et al., 2015), i.e. 115 mm 3 along R-L, A-P and S-I directions, ensuring realistic partial volume effects. Afterwards, we cropped the field-of-view along the foot-head direction to 200 mm (40 slices), in order to keep a tractable number of synthetic spinal cord voxels to analyse (i.e. 1700 voxels).
We used custom-written Matlab (The MathWorks, Inc., Natick, MA) code to synthesise signals for a rich multimodal quantitative MRI protocol encompassing of DW, qMT, inversion recovery (IR) and multi-echo time (multi-TE) imaging with shared imaging readout (protocol in Table 1, matching our rich in vivo MRI protocol). The total voxel-wise noise-free magnitude signal TOT was obtained as the weighted sum of the signals from WM, GM and CSF, i.e.  Above, is the 3  3 identity matrix, describes MT-weighting, = [0 0 1] T is aligned with the cord longitudinal axis and ( , T 1 , T 2 , AD, RD, , T 2 F , T 2 B , BPF) are tissuesspecific parameters, in this order: relative proton density (Mezer et al., 2013), macroscopic longitudinal and transverse relaxation rate (Smith et al., 2008), axial and radial diffusivity (Basser et al., 1994), free-to-bound pool exchange rate, free pool transverse relaxation rate, bound pool transverse relaxation rate, bound pool fraction (Henkelman et al., 1993). Eq. 5 models water relaxation as mono-exponential; diffusion as Gaussian, described by an axially symmetric diffusion tensor with primary diffusion direction aligned with the cord longitudinal axis; exchange between free and bound (i.e. myelin) protons according to the two-pool MT model (Henkelman et al., 1993). The MT-weighting factor was calculated via direct numerical integration of the two-pool Bloch equations (details in Supplementary Material S1), assuming a super-Lorentzian line shape for bound protons and simulating off-resonance pulse trains made of 25 sinc-Gaussian pulses (bandwidth: 122 Hz), each lasting 15 ms and with inter-pulse delay of 15 ms, as used before in spinal cord application (Battiston et al., 2018a).
We synthesised a unique noise-free signal profile in each tissue voxel by simulating within-tissue variability in WM and GM. This ensures that each synthetic voxel has its own unique sources of signal, avoiding obvious redundancies within the set of synthetic signals, as these could lead to overestimation of the performances of MP-PCA denoising (Ades-Aron et al., 2018). In practice, we drew voxel-wise values for each of ( , T 1 , T 2 , AD, RD, , T 2 F , T 2 B , BPF) from a tissue-specific Gaussian distribution, with parameters inspired by values known from literature (Battiston et al., 2018a;Grussu et al., 2015;Smith et al., 2008) (parameters in Table 2).

Denoising
We corrupted the synthetic signals with Gaussian noise at different signal-to-noise ratios (SNRs) (300 unique noise realisations on 1700 voxels), ranging from 10 to 40 (SNR evaluated with respect to the = 0 signal in WM for the DW measurements).
Afterwards, we used the Matlab implementation of the MP-PCA algorithm (http://github.com/NYU-DiffusionMRI/mppca_denoise) to denoise the synthetic spinal cord images at the various SNRs.
For our simulations, we processed MP-PCA matrices constructed by arranging spinal cord voxels within an individual MRI slice along rows and different MRI measurements along columns (i.e. slice-by-slice cord denoising). We implemented three different denoising strategies: 1. individual denoising of each modality among DW, IR, multi-TE and quantitative MT imaging respectively (importantly, IR, multi-TE have limited redundancy and would not theoretically qualify for MP-PCA, which is expected to remove little to no noise); 2. joint denoising of all modalities concatenated as one large set of measurements; 3. joint denoising of DW imaging concatenated with each of IR, multi-TE and qMT imaging in series respectively, which would be useful to describe cases when only one modality other than DW imaging is acquired.

Analysis
We evaluated the performance of MP-PCA denoising by studying the percentage relative error between the denoised signals TOT,denoised and the ground truth signal TOT . We estimated accuracy and precision of the different denoising strategies by calculating respectively the median of (such that the closer to zero, the higher the accuracy) and interquartile range (IQR) of (such that the lower, the higher the precision) within the synthetic spinal cord over the 100 noise instantiations. We also studied the distributions of normalised residuals ( ), evaluated according to Eq. 1.

In vivo study
We performed clinically viable, multi-contrast spinal cord qMRI scans on healthy volunteers and analysed them to characterise the performance of MP-PCA denoising on different qMRI modalities, devising denoising strategies for acquisitions with different levels of redundancy. The experimental sessions were approved by local research ethics committees.
Our qMRI protocols exhibit a unified signal readout, which is based on spin echo EPI, a typical choice for DW imaging. The shared readout ensures comparable noise characteristics across the different qMRI modalities, thus enabling joint denoising of different qMRI contrasts. The MRI protocol in vendor 1 encompasses DW, qMT, IR and multi-TE imaging, while in vendor 2 includes DW and multi-TE imaging. The MRI protocol in vendor 2 is less rich due to practical availability of pulse sequences. Nonetheless, it suffices to demonstrate the potential of joint multi-contrast denoising of modalities with different redundancies, and is representative of protocols required in multi-contrast techniques such as TEDDI .
In all systems, MRI scans were performed axially-oblique at the level of the cervical cord, with filed-of-view centred at the C2-C3 intervertebral disk (foot-head coverage of 60 mm).

MRI: vendor 1
The protocol developed on a 3T Philips Achieva machine, located at the UCL Queen Square Institute of Neurology (London, UK) consisted of multi-contrast, singleshot spin echo EPI scans with unified signal readout based on reduced field-of-view ZOOM technology (Wheeler-Kingshott et al., 2002), which enable 4 contrast mechanisms to be exploited: DW imaging, qMT imaging, IR imaging and multi-TE imaging (mTE, i.e. acquisitions of single-shot images at different TE). Salient sequence parameters, including information on b-values, echo/inversion times, off resonance saturation and cardiac gating are reported in Table 3.
The protocol also included an anatomical 3D FFE scan (flip angle of 7, TE of 4.1 ms, TR of 20 ms, resolution of 0.75  0.75  5 mm 3 and field-of-view of 180  240  60 mm 3 along R-L, A-P, S-I directions; ProSet fat suppression, 3 signal averages, scan time of 3 min : 30 s) and standard B0 and B1 field mapping for accurate qMT analysis. Both B0 and B1 mapping were based on 3D FFE acquisitions with resolution of 2.25  2.25  5 mm 3 and field-of-view of 215  206  60 mm 3 along R-L, A-P, S-I directions. B0 mapping was performed with the double-echo method (Jezzard and Balaban, 1995), with parameters: flip angle of 25, TE of 6.9 ms and 9.2 ms, TR of 50 ms, scan time of 1 min : 40 s. B1 mapping was instead performed via actual flip angle imaging (Yarnykh, 2007), with parameters: flip angle of 60, TE of 2.5 ms, TR of 30 ms, TR extension of 120 ms, scan time of 1 min : 40 s).
The nominal acquisition time was roughly 47 min, with variations depending on subject's heart rate. We scanned 4 healthy volunteers twice (2 males, age range 28-40), with the rescan performed within one month of the first scan.

MRI: vendor 2
For vendor 2, we performed scans on two separate 3T Siemens Prisma systems, located at the New York University School of Medicine (USA) and at the Neuroimaging Functional Unit of the University of Montreal (Canada).
The protocol consisted in exploiting 2 contrast mechanisms including DW imaging and multi-TE imaging with unified readout based on syngo ZOOMit reduced field-of-view technology (Rieseberg et al., 2002) (salient parameters including bvalues, TEs and cardiac gating are reported in Table 4). The protocol also included a 3D MEDIC scan for anatomical depiction (flip angle of 30, TE of 15 ms, TR of 625 ms, resolution of 0.50  0.50  5 mm 3 and field-of-view of 128  128  60 mm 3 along R-L, A-P, S-I directions; 3 signal averages, scan time of 6 min : 24 s).
The total scan time was 18 min : 57 s in the New York Prisma and 22 min : 19 s in the Montreal Prisma, with the scan time difference due to slightly higher number of diffusion directions being acquired in Montreal. Two subjects were scanned in New York (1 male, 28 years old; 1 female, 25 years old) and one subject (male, 28 years old) in Montreal after obtaining informed written consent. The vendor-provided 64 channel head-neck coil was used in both cases for signal reception.

Denoising
We implemented the same denoising strategies as in simulations: 1. individual denoising of each modality separately; 2. joint denoising of all modalities together; 3. joint denoising of DW imaging concatenated with each of IR, multi-TE and qMT imaging in series (multi-TE only for Prisma).
We performed denoising slice-by-slice to account for the anisotropic voxel-size and to limit the effect of potential between-shot signal fluctuations due to physiological noise (Summers et al., 2006). We proceeded as follows:  the spinal cord was identified on the mean DW image with SCT sct_propseg (De Leener et al., 2014);  all cord voxels of an MRI slice were arranged as one matrix and denoised with MP-PCA;  noise floor (Gudbjartsson and Patz, 1995) was subsequently mitigated on the denoised signals with the method of moments (Koay and Basser, 2006).

Post-processing
We performed motion correction on the concatenation of all acquired EPI images within an MRI session. Practically, we ran slice-wise rigid motion correction with sct_dmri_moco on the non-denoised scans, treating qMT, IR and multi-TE images as b = 0 scans. The estimated registration transformations were stored and used to correct all the denoised versions of each qMRI modality, as well as the nondenoised data. This was done to focus our analysis on the effect that thermal noise removal has on qMRI metrics.
Afterwards, we segmented the whole cord and the grey matter in the anatomical spinal cord scan respectively with sct_propseg and with sct_deepseg_gm. We also segmented the spinal cord in the mean DW EPI image with sct_propseg.
Lastly, we co-registered the anatomical spinal cord scan to the mean EPI image with sct_register_multimodal, using dilated spinal cord masks in the two image spaces to guide registration (dilation performed with NifTK seg_maths, available at http://github.com/NifTK/NifTK ). The estimated warping field transformation was used to warp the grey matter mask to EPI space, which was subsequently used to obtain a white matter mask by subtracting it from the whole-cord mask. For vendor 1, the warping field was also used to resample the B0 and B1 magnetic field maps to the EPI space for downstream model fitting.

Analysis
We calculated the normalised signal residuals (Eq. 1) for all denoising approaches and for all subjects, sessions and vendors, and evaluated their distributions within the spinal cord.
We also characterised values of all qMRI metrics by calculating the median within grey and white matter for all denoising strategies (including no denoising).
Finally, we quantified each metric variability by calculating a percentage coefficient of variation (CoV) within grey matter and within white matter for all denoising strategies (including no denoising). We defined CoV as where iqr is the interquartile range of a metric within grey/white matter, measuring the metric variability, while is the median value of the metric within the same tissue. We hypothesise that effective denoising would reduce noise-induced metric variability, resulting in lower iqr and unchanged median and hence lower CoV, under assumption that variability due to noise is much larger than the biological variability (please see figure 4 of (Ades-Aron et al., 2018)). Fig. 1 shows distributions of residuals (in the semilog scale, plotted against the squared residuals) from simulations. These are reported for the different qMRI modalities and for different denoising strategies. In such plots, Gaussian residuals align along a straight line. The charts reveal that residuals are Gaussian for at least 3 . Moreover, they also highlight that joint multi-contrast MP-PCA denoising mitigates more efficiently noise for those modalities featuring limited redundancy than if denoised independently, such as IR imaging and mTE imaging: their residuals are closer to the ideal unit Gaussian when denoised jointly with more redundant modalities than when denoised alone. In contrast, denoising performance does not improve for joint multi-contrast MP-PCA of those modalities that are intrinsically highly redundant (i.e. DW and qMT imaging). Fig.2 shows percentage relative error accuracy (top row: error median) and precision (bottom row: error IQR) of the denoised signals compared to the noise-free ground truth, for different qMRI modalities and different denoising strategies. While plots do not highlight any noticeable differences in terms of accuracy for the different denoising strategies (i.e. joint denoising or individual denoising, since their confidence intervals overlap perfectly), they do suggest that better precision (i.e. error IQR closer to zero) can be achieved for modalities that are intrinsically limited in redundancy, when these are denoised jointly with more redundant modalities. For example, IQR drops from 8% to 5% at SNR = 10 for mTE when it is denoised jointly with all modalities, as compared to when mTE is denoised along. Similar to what is reported in Fig. 1, no appreciable improvement of denoising performance is observed with joint multi-contrast denoising for modalities that intrinsically feature high redundancy. This is apparent for qMT and even more so for DW imaging, since their percentage error IQR does not change when these are denoised jointly with other modalites).  Fig. 7; DW and mTE imaging) for different denoising strategies. Visual inspection suggests that MP-PCA denoising generates less noisy maps, especially for vendor 1. The most striking examples of improved parameter estimation are seen in both vendors for DW imaging parameter MK. Additionally, improvements on visual inspection are apparent for other qMRI metrics such as BPF and T 1 , especially for joint multimodal denoising.

In vivo study
Tables 5 and 6 report median values of all qMRI metrics in grey and white matter for the different denoising strategies (Table 5: vendor 1; Table 6: vendor 2, pooling together results from the two systems). The tables reveal contrasts between grey and white matter in various metrics. Examples that are consistent between vendors include: higher FA and MD in white compared to grey matter; similar MK in grey/white matter; slightly higher T 2 in white compared to grey matter. Other examples from vendor 1 include: similar BPF and T 1 in grey/white matter; higher exchange rate in grey compared to white matter. The tables also show that systematic differences between the data sets acquired with the two vendors exist, as for example: higher T 2 and MK and lower MD in data from vendor 2 compared to 1; different grey/white matter contrasts in FA. Notably, Tables 5 and 6 also demonstrate that denoising introduce little to no biases in the quantitative parametric maps. In all cases and for both vendors the tissue-wise medians never differ for more than 5% compared to the values obtained without any denoising.
Tables 7 and 8 report within-grey and within-white matter CoV for the various qMRI metrics and for different denoising strategies. Table 7 reports figures from vendor 1,  while Table 8 from vendor 2 (data from both systems from vendor 2 pooled together). The tables show that MP-PCA denoising leads to reductions of CoV for various metrics of 5% or more compared to the case with no denoising, as for example for FA, MK, BPF and T 1 for vendor 1 and MK for vendor 2. Some increases of CoV are observed (for example for MD in white matter for vendor 2). For vendor 1, the strongest reductions in CoV are observed for joint multimodal MP-PCA denoising.

Summary and key findings
This work demonstrates the advantages of multimodal qMRI of the spinal cord in vivo with unified MRI signal readout. The unified readout enables matching resolution and distortions across different contrasts, thereby also facilitating joint analyses and computational modelling of multi-contrast signal. These include MP-PCA denoising, a recent noise mitigation technique tailored for redundant qMRI protocols, i.e. such that the number of measurements is much higher than the underlying number of signal components. The unified readout enables efficient MP-PCA denoising of modalities that would exhibit little redundancy if denoised individually.
Our key findings are that a unified readout enables reliable and detailed microstructural characterisation of the human cervical spinal cord in clinical setting, providing metrics of relaxometry and diffusion as well as myelin-sensitive indices with matched resolution and distortions. Moreover, MP-PCA appears as a valid tool to improve the intrinsic quality of unified readout acquisitions, as supported by both in vivo and in silico data. Finally, this approach is feasible on 3T MRI systems from two major vendors.

In silico study
We have designed and run computer simulations to test whether a unified readout offers opportunities for MP-PCA denoising of qMRI modalities that exhibit limited redundancy (i.e. a number of measurements comparable to the number of signal generators), for which effective MP-PCA denoising remains challenging. For this purpose, we synthesised spinal cord qMRI data for a rich protocol consisting of DW, qMT, IR and mTE imaging, using anatomical information from the SCT atlas. We generated signals using literature values for microstructural parameters in MRI signal forward models, ensuring within-tissue variability to avoid obvious signal redundancy.
Our simulations suggest that a unified readout has indeed the potential of supporting more efficient MP-PCA denoising for modalities limited in redundancy, as for example mono-exponential and multi-exponential (Does et al., 2019) relaxation mapping. Denoising these modalities jointly with more redundant modalities enables more efficient noise mitigation in the former, thus suggesting that some of the information carried by the signal generators may be at least in part shared among MRI contrasts. Interestingly, joint multimodal denoising did not affect those modalities that are already redundant, as for example DW imaging, thus providing another piece of evidence supporting the idea that different contrasts may share some of the signal sources.

In vivo study
In this paper we tested multi-parametric qMRI of the spinal cord with unified readout using 3T MRI scanners from two major vendors (Philips and Siemens). We obtained a number of qMRI metrics that are promising biomarkers in several spinal cord conditions, including diffusion parameters FA, MD and MK; qMT two-pool BPF and exchange rate ; relaxation time constants T 1 and T 2 . We studied whether MP-PCA improves the quality of raw MRI signals as well as quantitative maps, by evaluating distributions of signal residuals as well as within-tissue variability of metrics via a CoV.
Our multi-vendor data demonstrate the feasibility of implementing reliable multiparametric qMRI of the spinal cord with unified readout. A unified readout provides matched resolution and distortions across MRI contrasts, and ensures comparability of signals across a rich set of qMRI measurements. Moreover, it enables the development of unified analysis pipelines, spanning from motion correction, to data denoising and potentially model fitting, paving the way to joint modelling of multicontrast signals (Kim et al., 2017). Importantly, it may be useful in techniques that combine information from diffusion with relaxation/myelin-sensitive indices, as for example g-ratio MRI (Campbell et al., 2018;Duval et al., 2017;Stikov et al., 2015)), where matched EPI distortions (Irfanoglu et al., 2015) are crucial (Campbell et al., 2018). Our unified readout provides images with the same distortions across all contrasts and without significant losses of resolution for scans that would be typically performed with a different readout compared to DW imaging (Duval et al., 2017;Ljungberg et al., 2017), i.e. myelin mapping and relaxometry.
In this paper we also exploited synergies across MRI contrasts to evaluate the performance of MP-PCA denoising in the spinal cord. To our knowledge, this is the first time that the practical utilities of MP-PCA has been characterised in detail in the human spinal cord in vivo. Moreover, this is the first time that MP-PCA has been tested across that many MRI contrasts in a unified acquisition. Our analyses demonstrate that MP-PCA effectively mitigates noise in all modalities and for both vendors. Importantly, distributions of residuals show that the efficiency of MP-PCA in enhancing the quality of modalities with limited redundancy (i.e. IR and mTE imaging) can be improved by denoising these modalities jointly with more redundant schemes. Such results mirror findings in simulations, and thus suggest that some of the information about MRI signal sources may be shared across contrasts.
Our joint multimodal denoising relies on the hypothesis of noise homoscedasticity across MRI contrasts. Supplementary material S2 and S3 shows that the estimated noise level on modalities other than DWI follows the same trends as those of estimates from DWI in both simulations (S2) and in vivo (S3). Nonetheless, the two supplementary documents also demonstrate that estimating the noise level is a very challenging task: the estimates of the noise level are highly variable per se. Moreover, Supplementary material S3 also reveals systematic differences between noise standard deviation estimates from DW imaging compared to other modalities, such as qMT. This is likely attributed to the stronger departures from the hypothesis of Gaussian noise underlying MP-PCA in DW imaging, due to lower SNR and stronger noise floor effects (Koay and Basser, 2006), and to the fact that qMT suffers from stronger physiological noise that may resemble thermal noise (qMT is not cardiac gated). Nonetheless, it should be remembered that noise levels estimated on modalities with limited redundancy (e.g. multi-TE imaging) are not reliable, as the limited number of measurements does not allow the MP distribution to emerge, spoiling the detection of noisy eigenvalues in MP-PCA (Veraart et al., 2016b). Importantly, such differences in terms of noise level estimates among modalities introduce little to no bias in downstream quantitative parameter maps, and therefore do not appear to be a concerning issue for practical MP-PCA deployment.
We also investigated the effect of MP-PCA denoising on the quality of popular parametric maps. To this end, we studied median values of metrics within grey/white matter as well as metric variabilities as quantified by a CoV. Our experiments show that MP-PCA introduces little to no biases in any of the metrics, irrespectively of the chosen denoising strategy (joint multimodal denoising vs modality-wise denoising). The difference in median values between metrics obtained with denoising compared to the case with no denoising are 5 % or less. These differences, which are very low, likely reflect the intrinsic susceptibility of the different model fitting routines to noise fluctuations and noise floors, and are therefore expected since noise-floor mitigation (Koay and Basser, 2006) was performed following MP-PCA denoising. Conversely, MP-PCA does decrease metric variability, at it leads to tissue-wise CoV reduction much higher 5% for many metrics, as for example for MK (both vendors), FA, BPF and T 1 (vendor 1). The reduction in variability is the highest for metrics like MK, which carry important information about tissue microstructure, and that are notoriously difficult to estimate (Veraart et al., 2011).
Our parametric maps follow known trends and contrasts, with some differences in terms of relaxometry metrics, e.g. low contrast white/grey matter contrast for T1 and T2. This difference may be explained by residual CSF pulsation that corrupts neighbouring white matter signals, and by the fact that literature values for T1 and T2 are typically obtained with different readout strategies compared to the employed single-shot EPI (Smith et al., 2008). Another explanation, especially for vendor 1, may be related to the coarse resolution of the anatomical scan, required to limit scan time, as this was used for grey matter segmentation potentially introducing partial volume effects in the tissue masks. Overall, while grey/white matter contrasts in parametric maps are similar in data from both vendors, systematic differences between metric values (Table 5 vs 6) and variability (Table 7 vs 8) are seen. Several factors may have contributed to these differences between vendors, namely in: intrinsic SNRs; rFOV techniques; resolution of the anatomical scan used for grey/white matter segmentation; parallel imaging/reconstruction techniques; qMRI protocol; betweensubject biological differences.
Finally, we point out that we took care to use the same registration transformations to correct for motion in all denoising strategies, estimating motion on the non-denoised data. We followed this motion correction strategy on purpose, as our focus was to study the effects of thermal noise mitigation on parametric maps. It is possible that the benefits of MP-PCA may extend beyond thermal noise mitigation and may also improve post-processing such as motion correction, as shown in other studies (Ades-Aron et al., 2018), which will be the subject of future investigations.

Considerations and future directions
We acknowledge a number of limitations of our approach.
Firstly, our unified protocol was more comprehensive on the system from vendor 1 as compared to vendor 2. This was due to the practical availability of MRI sequences at the time of acquisition of the data. In future we plan to expand our protocol on vendor 2 to better characterise MP-PCA performance across contrasts.
Secondly, the DW imaging protocol for the vendor 2 differed between the system in New York and the one in Montreal, with the latter being slightly longer. This was due to a choice in the design of the protocol in Montreal, which would enable the inclusion of the scan in other ongoing group studies. Nonetheless, the acquisition suffices to demonstrate the potential of a unified readout and to explore multimodal denoising.
Furthermore, in the future we plan to expand our sample size to better characterise the potential of our unified acquisition and denoising in real clinical settings on patients as, for example, in multiple sclerosis.

Conclusions
Multi-parametric qMRI of the spinal cord with unified readout (i.e. with matched resolution and distortions) is advantageous and provides robust microstructural metrics sensitive to axonal characteristics, relaxometry and myelin. Our unified acquisition paves the way to joint modelling of multi-contrast signals, and offers unique opportunities for image quality enhancement with techniques such as MP-PCA denoising, proven here to be a useful pre-processing step in spinal cord MRI analysis pipelines.

Data and code availability statement
The synthetic spinal cord phantom is made freely available online (http://github.com/fragrussu/PaperScripts/tree/master/sc_unirea dout/sc_phantom). All scripts written to analyse the data are also made openly available (http://github.com/fragrussu/PaperScripts/tree/master/sc_unirea dout). The in vivo data cannot be made openly available online due to privacy issues of clinical data according to GDPR regulations. Researchers interested in accessing the in vivo data from vendor 1 can contact Prof Claudia Gandini Wheeler-Kingshott (c.wheeler-kingshott@ucl.ac.uk). A data sharing agreement enabling academic and research use will be stipulated. Researchers interested in accessing the in vivo data from vendor 2 can contact: Prof Timothy Shepherd (timothy.shepherd@nyulangone.org) for the New York data; Prof Julien Cohen-Adad (jcohen@polymtl.ca) for the Montreal data. A data sharing agreement will be stipulated with either New York University, Polytechnique Montreal or both, enabling academic and research use.
The code for qMT signal synthesis and analysis (two-pool fitting) is not openly available online. Researchers interested in accessing the code can contact Dr Marco Battiston (marco.battiston@ucl.ac.uk), who would release a copy under a license/sharing agreement enabling academic and research use.  Table 1 Sequence parameters used to simulate synthetic multimodal spinal cord scans. In the table, DW, qMT, IR and multi-TE stand respectively for diffusion-weighted, quantitative magnetisation transfer, inversion recovery and multi-echo time. All of DW, qMT, IR and multi-TE imaging rely on the same spin echo EPI readout with long TR (i.e. such that it is hypothesised that TR >> T1). For qMT, each of the 4 repetitions of 11 MTweighted measurements is characterised by a different delay between the end of the off-resonance pulse train and the readout, i.e. {17, 95, 173, 251} ms. The offresonance pulse train in qMT was made of 25 sinc-Gaussian pulses (bandwidth: 122 Hz), each lasting 15 ms and with inter-pulse delay of 15 ms (Battiston et al., 2018a Table 2 Tissue parameters used to generate the synthetic spinal cord scans. Values are inspired by previous literature (Battiston et al., 2018a;Grussu et al., 2015;Smith et al., 2008). For white/grey matter, within-tissue variability was simulated by drawing parameter values from a Gaussian distribution and assigning the obtained values to different voxels. The mean and standard deviation of the Gaussian distributions are reported in the table (standard deviation within brackets, equal to 10% of the mean). For the cerebrospinal fluid (CSF), tissue parameters were fixed to the same values across all CSF-containing voxels.  Table 3 Salient sequence parameters for the qMRI protocol with unified readout implemented on vendor 1 (Philips Achieva, London, UK). DW, qMT, IR and multi-TE stand for diffusion-weighted, quantitative magnetisation transfer, inversion recovery and multiecho time. Consistently with simulations, in qMT each repetition is characterised by a different delay between the end of the off-resonance train and the readout, i.e. {17, 95, 173, 251} ms. qMT off-resonance trains were made of 25 sinc-Gaussian pulses (bandwidth: 122 Hz), each lasting 15 ms and with inter-pulse delay of 15 ms (Battiston et al., 2018).     Table 7 Percentage CoV in grey and white matter for the various qMRI metrics obtained from vendor 1 (London, UK) with different denoising strategies. The table reports CoV = 100 iqr/median, where iqr and median are respectively the interquartile range and the median of a metric within grey/white matter. CoV from different subjects and scans are pooled so that figures report mean and standard deviation (in brackets) of CoV across subjects and scans. Reductions in mean CoV with more than 5% as compared to the case with no denoising are labelled in bold font and grey shadowing (no increase of mean CoV greater than 5% is observed).  Table 8 Percentage CoV in grey and white matter for the various qMRI metrics obtained from vendor 2 (New York, NY, USA and Montreal, Canada) with different denoising strategies. The table reports CoV = 100 iqr/median, where iqr and median are respectively the interquartile range and the median of a metric within grey/white matter. CoV from different subjects and scans are pooled so that figures report mean and standard deviation (in brackets) of CoV across subjects and scans. Improvements of mean CoV greater than 5% compared to the case with no denoising (i.e. lower CoV) are labelled in bold font and grey shadowing (no worsening of mean CoV greater than 5% is observed).  Fig. 1. Distributions of normalised residuals from simulations for different denoising strategies. Each plot reports the logarithm of the residual histogram, scattered against the corresponding values of squared residuals. In these plots, Gaussian residuals appear as straight lines. First row (panels A to D): simulated SNR of 10; second row (panels E to H): simulated SNR of 15; third row (panels I to L): simulated SNR of 20; fourth row (panels M to P): simulated SNR of 40. SNR is evaluated with respect to the synthetic white matter signal in the diffusion scan at b = 0. Different columns refer to different MRI scans: from left to right, DW imaging (panels A, E, I, M), qMT imaging (panels B, F, J, N), IR imaging (panels C, G, K, O), mTE imaging (panels D, H, L, P).

Fig. 2.
Accuracy and precision of different denoising strategies as obtained from percentage relative errors (percentage errors between denoised signals and noisefree ground truth signals) in simulations. Panels A to D (top) show median percentage at different SNR levels, and represent a measure of accuracy (the closer to zero, the higher the accuracy; DW imaging in A, qMT imaging in B, IR imaging in C, mTE imaging in D). Panels E to H (bottom) show percentage relative error interquartile ranges at different SNR levels, and represent a measure of precision (the closer to zero, the higher the precision; DW imaging in E, qMT imaging in F, IR imaging in G, mTE imaging in H).

Fig. 3.
Examples of MP-PCA denoising in one subjects who was scanned with vendor 1. Panels A, C, E, G show raw and denoised images, obtained according to different strategies; panels B, D, F, H show the corresponding distributions of normalised residuals in log-scale, plotted against the squared residuals (in these plots, Gaussian residuals align along a straight line). DW imaging: panels A (images) and B (residuals); qMT imaging: panels C (images) and D (residuals); IR imaging: panels E (images) and F (residuals); mTE imaging: panels G (images) and H (residuals). Ant, Post, Right, Left respectively indicate subject's anterior, posterior parts and right and left sides.  Different rows illustrate the metrics obtained according to different denoising strategies (no denoising; independent denoising of each modality; various combinations of joint multi-modal denoising). Quantitative maps are overlaid onto the mean non-DW image and shown within the cord only. The same anatomical conventions with regard to subject's anterior, posterior parts and right and left sides as in Fig. 3 are used.  Fig. 6. Examples of quantitative maps from vendor 2 (Siemens Prisma system located in New York, USA). From top to bottom: FA, MD, MK (DW imaging); T2 (mTE). Different rows illustrate the metrics obtained according to different denoising strategies (no denoising; independent denoising of each modality; various combinations of joint multi-modal denoising). Quantitative maps are overlaid onto the mean non-DW image and shown within the cord only. The same anatomical conventions with regard to subject's anterior, posterior parts and right and left sides as in Fig. 4 are used.  Fig. 7. Examples of quantitative maps from vendor 2 (Siemens Prisma system located in Montreal, Canada). From top to bottom: FA, MD, MK (DWI); T2 (mTE). Different rows illustrate the metrics obtained according to different denoising strategies (no denoising; independent denoising of each modality; various combinations of joint multi-modal denoising). Quantitative maps are overlaid onto the mean non-DW image and shown within the cord only. The same anatomical conventions with regard to subject's anterior, posterior parts and right and left sides as in Fig. 4 are used.

Supplementary material S1: the two-pool MT model
The two-pool magnetisation transfer (MT) model is used in this paper in both simulations and in vivo analyses. The model describes the evolution of the magnetisation of two exchanging 1 H pools, i.e. free protons in bulk water and protons bound to macromolecules (Henkelman et al., 1993), in presence of off-resonance irradiation 1,off ( ) with carrier frequency 0 +  (with  being the offset with respect to water resonance frequency 0 ), assumed to be played out at times 0 ≤ ≤ off . On-resonance excitation is then assumed to be played out at > off , after a generic delay delay . Let T the free and bound proton magnetisations during the offresonance irradiation, and with T 1 F and T 1 B and with T 2 F and T 2 B their respective longitudinal and transverse relaxation times (Portnoy and Stanisz, 2007). Let us also indicate with the exchange rate between free to bound protons, with BPF the bound pool fraction (fraction of protons belonging to the bound pool) and with the proton gyromagnetic ratio ( 2 ≅ 42.577 MHz T ).
The MT-weighting factor in Eq. 5 (main manuscript) is defined as the value of the free pool longitudinal magnetisation at the end of off-resonance irradiation, i.e.: = F ( = off ).
( 1.1) The value of in Eq. S1.1 is obtained directly by numerical integration of the two-pool Bloch equations, which are written as: ( 1.2) For the integration of equation S1.2, it is assumed that F ( = 0) ≝ [0 0 (1 − BPF)] T , B ( = 0) ≝ [ 0 0 BPF ] T , B ( = off ) = B ( = off ) ≈ 0 and that the bound pool off-resonance absorption can be modelled based on a super-Lorentzian line shape. Under this assumption, the absorption term RF B is a function of ( ,  , T 2 B ) and is modelled as (Portnoy and Stanisz, 2007): For the practical integration of Eq. S1.2, it is assumed that T 1 B ≈ 1 s (Battiston et al., 2018a), and that the free pool longitudinal relaxation time (i.e. T 1 F ) is linked to the observable T 1 relaxation time as (Portnoy and Stanisz, 2007): .

Supplementary material S3: estimated noise levels in vivo
In this supplementary material we compare the estimates of noise levels provided by MP-PCA denoising on different modalities (in vivo data).
In the following supplementary figures, we scatter the estimated noise standard deviation obtained from denoising the DWI against the values obtained by denoising the other modalities, either alone or jointly with DWI. Moreover, we report the scatter plots with and without noise floor mitigation, and for the two vendors.
Results are shown in both cases when denoised images and estimated noise level are corrected for Rician noise bias with the method of moments [1] and whey they are not. Fig. S3.1. Estimated noise standard deviation for vendor 1 (Philips Achieva) obtained denoising each modality (qMT, IR and mTE) independently. Values are scattered against the noise level estimated from the DWI alone. In panel A) (left), the estimated noise level has been corrected for Rician bias, while in panel B) (to the right) this correction has not been performed. When generating the plots, data from all scans of all subjects have been pooled together. Figure S3.2. Estimated noise standard deviation for vendor 1 (Philips Achieva) obtained denoising each modality (qMT, IR and mTE) jointly with DWI. Values are scattered against the noise level estimated from the DWI alone. In panel A) (left), the estimated noise level has been corrected for Rician bias, while in panel B) (to the right) this correction has not been performed. When generating the plots, data from all scans of all subjects have been pooled together. Fig. S3.3. Estimated noise standard deviation for vendor 2 (two Siemens Prisma located in New York, USA and Montreal, Canada) obtained denoising mTE independently from DWI. Values are scattered against the noise level estimated from the DWI alone. In panel A) (left), the estimated noise level has been corrected for Rician bias, while in panel B) (to the right) this correction has not been performed. When generating the plots, data from all scans of all systems have been pooled together. Fig. S3.4. Estimated noise standard deviation for vendor 2 (two Siemens Prisma located in New York, USA and Montreal, Canada) obtained denoising mTE jointly with DWI. Values are scattered against the noise level estimated from the DWI alone. In panel A) (left), the estimated noise level has been corrected for Rician bias, while in panel B) (to the right) this correction has not been performed. When generating the plots, data from all scans of all systems have been pooled together.