Reliability of multi-parameter mapping (MPM) in the cervical cord: A multi-center multi-vendor quantitative MRI study

,


Introduction
Quantitative MRI (qMRI) can measure microstructural features of the brain and spinal cord tissue. QMRI provides a mean to investigate tissue properties with higher sensitivity both in the healthy and diseased central nervous system ( Edwards et al., 2018 a;Leutritz et al., 2020 ;Weiskopf et al., 2021 ;Cohen -Adad et al., 2022 ). In theory, qMRI readouts are insensitive to differences in MR software and hardware and thereby increase the comparability across different time points and sites ( Weiskopf et al., 2021 ;Edwards et al., 2018 b;Voelker et al., 2021 ). This feature renders qMRI attractive to clinical trials that are most often run as multi-center studies. Thus, a pressing need exists to provide a standardized acquisition protocol of qMRI and image processing tools for the brain but also the spinal cord. www.nisci-2020.eu ) investigating a therapeutic effect based on the Anti-Nogo-A antibody treatment ( Freund et al., 2006 ) in acute SCI patients.
In this study, we aimed to additionally determine the intra-site repeatability (scan-rescan) and inter-site reproducibility of the cervical MPM data obtained in the same traveling heads study for which we previously reported the intra-site repeatability (scan-rescan) and inter-site reproducibility of the brain MPM data ( Leutritz et al., 2020 ).

Ethics statement
Informed written consent was obtained from each subject prior to each measurement per site, and all sites obtained local ethics approval from their institutional ethical review.

Participants and MRI sites
The study was conducted on six 3T MRI scanners with different hardware and software ( Table 1 ), having four Siemens MRI scanners (Erlangen, Germany) and two Philips MRI scanners (Eindhoven, Netherlands). Five healthy participants (2 females, age: 32.4 ± 6.0 years [mean ± standard deviation]) underwent a scan-rescan MRI session across all sites with an average interscan break of approximately two hours between the two measurements ( Leutritz et al., 2020 ). All study participants have been scanned using 3T MRI scanners equipped with head and neck RF coil.

MRI acquisition
The MPM protocol was implemented based on product sequences available on the manufacturer's clinical MRI systems. Threedimensional (3D) data acquisition comprised of three multi-echo spoiled gradient echo scans (i.e. fast low angle shot [FLASH] sequences on Siemens scanners and multi-echo fast field echo [mFFE] sequences on Philips scanners) with MT, T1, and PD contrast weighting. Additional reference scans for radio-frequency (RF) transmit (B1 + ) and receive (B1 − ) fields for bias correction within the hMRI-toolbox ( Tabelow et al., 2019 ) were acquired. The total acquisition time was 18:45 min on the Siemens scanner and 23:58 min on the Philips scanner. Therefore, we adapted a previous protocol used in an SCI longitudinal study ( Freund et al., 2013 ) by reducing the total number of acquired echoes, repetition time (TR) and maximum echo time (TE), which still allows for mapping R2 * well in subcortical -short T2 * -areas. Optimal flip angles were derived from the Ernst angle (based on respective TR and an average longitudinal relaxation time (T1) of 1050 ms for brain tissue) by multiplying them with a factor of 0.4142 for PD-and MT-weighted acquisitions, and by 2.4142 for T1-weighted acquisitions, respectively, in order to minimize noise propagation . The acquisition protocols shared the following common parameters across all platforms ( Table 2 ): TR of PD-and T1-weighted contrasts: 18 ms; flip angles for MT-, PD-and T1-weighted contrasts: 6°, 4°, 25°, respectively; six equidistant echoes; 1 mm 3 isotropic reconstruction voxel size; readout (RO) field of view (FoV): 256 mm; base resolution: 256 pixels in the RO direction; 176 slices; readout in the head-foot direction, inner phase encoding loop in the left-right ( "slice ", fast phase encoding) direction, outer phase encoding loop in the anterior-posterior direction ( "phase ", slow phase encoding); RO bandwidth: 480 Hz/pixel; elliptical k-space coverage; parallel imaging speedup factor of 2 in the slow/slice phase encoding direction. We accounted for differing RF spoiling characteristics between vendors' sequence implementations ( Table 2 ) in the processing. Since the MT pulse details cannot be changed on different vendors without pulse sequence programming and appropriate research agreement, which is not practical for the clinical scanners, we used standard settings ( Table 2 ). Additionally, a post-acquisition harmonization between vendor implementations was performed ( Leutritz et al., 2020 ). The minimum TR of the MT-weighted sequences was restricted at Philips scanners due to specific absorption rate (SAR) constraints and thus extended the total acquisition time in comparison to Siemens scanners. The B1 + field mapping methods differed across vendors and sites. At Siemens sites, vendor-supplied sequences were used. It was either based on spin-echo and stimulated echo acquisitions and is similar to the customized sequence by Lutti et al. (2010) or based on the gradient echo sequence with ultrafast turbo-FLASH readout (available from version VE11 onwards on Siemens scanners). At the Philips scanners a vendor based actual flip angle imaging (AFI) technique was used ( Yarnykh, 2007( Yarnykh, , 2010. Moreover, to obtain B1 − maps and correct for the apparent sensitivity changes between MPM acquisitions induced by potential head motions, we used low-resolution 3D spoiled gradient echo volumes once with RF head and neck coils and once with the body coils consequently prior to each MPM modality acquisition. The ratio between RF coil sensitivity maps results into a relative net RF received field sensitivity (B1 − ) ( Papp et al., 2016 ;Tabelow et al., 2019 ). For more details in the B1 + and B1 − mapping acquisition parameters please refer to the previous published report based on the brain MPM ( Leutritz et al., 2020 ).

Data quality control
The acquisition parameters of each scan (as stored in the DICOM header) were manually checked against standard settings to detect inconsistencies in the data acquisition ( Leutritz et al., 2020 ). Intermediate data volumes, segmentations, and parameter maps were systematically checked visually, to detect misregistration or wrong scaling of the quantitative maps.

Creating MPM maps and cross-sectional cord area
The MPM data were processed using the hMRI-toolbox ( www.hmri.info ,versions v0.1.1-beta and v0.1.2-beta) ( Tabelow et al., 2019 ) within the SPM12 framework (version 6906; FIL, London, UK) in MATLAB (version R2017b; Mathworks, Natick, Massachusetts, USA on GNU/Linux computers with x86_64 architecture) ( Fig. 1 A). The main processing steps included the data conversion, calculation of quantitative maps and reproducibility analyses. For optimal segmentation and registration of volumes, we first applied auto-realignment as implemented in the hMRI-toolbox ( Tabelow et al., 2019 ). R1, PD, R2 * and MT maps were estimated from the multi-echo data in combination with the B1 + and B1 − measurements using the "Create hMRI maps module " in the hMRI toolbox and applying the probabilistic atlas of the brain and neck ( Blaiotta et al., 2018 ;Tabelow et al., 2019 ;. Estimation of R1, PD and MT maps was based on the rational approximation of the signal  with additional corrections for B1 + inhomogeneities ( Tabelow et al., 2019 ). The R2 * maps were calculated based on all echoes acquired across all contrasts using the ESTATICS estimation scheme ( Weiskopf et al., 2014 ). To correct for R2 * -and remaining B1 − -related bias of the PD maps, the signal was extrapolated to TE = 0 ( Leutritz et al., 2020 ).
To extract the quantitative maps of the cervical cord we used the spinal cord toolbox (SCT) ( De Leener et al., 2017 ). The spinal cord from levels C1 to C3 was automatically segmented based on the MT maps. Next, MT maps were co-registered with the MNI-Poly-AMU template ( De Leener et al., 2018 ). Using the spinal cord segmentations and vertebral labels mean values of MT, PD, R1, and R2 * were extracted from the parameter maps for each vertebral level (C1-C3) separately.
Spinal cord cross-sectional area: To estimate the spinal cord crosssectional area, the T1-weighted images of the MPM were segmented us-ing the SCT ( De Leener et al., 2017 ). First, an averaged T1-weighted MRI was calculated across six echoes to increase the contrast between the spinal cord and CSF signal. The spinal cord was then automatically segmented on the averaged images using sct_deepseg_sc ( Gros et al., 2019 ) and co-registered with the MNI-Poly-AMU ( De Leener et al., 2018 ). Finally, the cross-sectional area (CSA), left-right width (LRW), and anterior-posterior width (APW) of the spinal cord were estimated based on the spinal cord segmentation for each vertebral level (C1-C3) ( Fig. 1 B).

Analysis of inter-and intra-site reproducibility
To determine intra-and inter-site reproducibility of the MPM metrics and the macrostructural parameters (CSA, LRW, APW) of the spinal cord, coefficients of variance (CoV) within and between sites were calculated across C1 to C3 levels -for each MPM map and each macrostructural parameter ( Fig. 1 ). To assess systematic bias, mean parameter values were additionally compared between sites. The intra-site CoV was determined as the standard deviation ( intra ) of the extracted parameter (average over all voxels within each ROI) estimated from scan and rescan data over the mean (μ intra ) of both parameters at a single site for each vertebral level separately: The inter-site CoV was determined as the standard deviation ( inter ) over the mean (μ inter ) across all scans for a specific subject, comprising bias observed at individual sites: ROI is C1, C2 and C3 level of the cervical cord, site is the scanning center where the data were acquired, subj is the study participant identifier, and parameter is the quantitative parameters such as MT, R1, R2 * , PD, CSA, LRW, and APW. This represents the precision of the metric within the same subject.

Inter-site and vendor-bias
The site-specific bias ( Δ) was also defined as the normalized difference of the mean (μ intra ) of both maps at a single site minus the mean of all μ intra across all sites (normalized by the μ intra across all sites). For

Fig. 2.
Mean of MPM maps across all sites across levels C1-C3 of cervical cord. The mean per subject was calculated across the scan and rescan measurements as well as across sites. Box plot = interquartile range (q3 -q1); bold horizontal line = median, feathers = data range of values across subjects. Cross = outlier, which exceeds the range of (q1-1.5 × [q3 − q1]; q3 + 1.5 × [q3 − q1]), where q1 and q3 are the 25th and 75th percentiles of the sample data, respectively. MT: Magnetization transfer saturation, PD: proton density, R1 and R2 * : longitudinal and effective transverse relaxation, respectively. better assessment of a systemic bias caused by different MRI systems, the site-specific bias was averaged across ROIs (C1-C3): To assess the differences due to slightly different acquisition MPM protocols across different vendors as well as hardware and software available at each site, (e.g., different RF pulses or spoiling characteristics), the data were sub-grouped as follows: (a) data from all sites equipped with Siemens scanners, (b) data from all sites equipped with Philips MRI scanner, (c) data from all sites, irrespective of manufacturer. For each of these subgroups, the intra-site CoV, inter-site CoV, and site-specific bias were determined according to the same methods described above. Next, the root-mean-square (RMS) value (x RMS ) was calculated for each of the measures (intra-site CoV, inter-site CoV and the site-specific bias) to average the different subjects, sites, and ROIs, respectively, using Eq. (4 ). Thereby, a single output measure was obtained for each subgroup, which allows for better comparison between groups.

Mean values of the MPM maps
The mean values of the MT, PD, R1 and R2 * maps of the cervical cord across vertebral levels (C1-C3) are shown in Fig. 2 . The mean value of MT across C1-C3 cervical cord was 3.28 ± 0.45 p.u. (percent unit), the mean value of PD was 71.98 ± 4.94 p.u. (percent unit), the mean value of R1 was 0.88 ± 0.09 s − 1 , and R2 * mean was 19.98 ± 3.00 s − 1 . The value of the CSA based on the T1-w MRI was 62.04 ± 4.85 mm 2 , the APW and LRW were 7.28 ± 0.34 mm and 10.79 ± 0.84 mm, respectively ( Fig. 3 ).

Intra-and inter-site reproducibility
The mean intra-site CoV across cervical levels C1 to C3 was 7.3% for MT, 2.5 % for PD, 5.4% for R1 and 7.5% for R2 * across the C1-C3 levels based on ROI analysis ( Fig. 4 A-D). The mean inter-site CoV across cervical levels C1 to C3 was 10.5% for MT, 6.1% for PD, 7.7% for R1 and 12% for R2 * ( Fig. 5 ). The mean intra-site CoV of levels C1 to C3 of the cervical cord of CSA was 2.0% and 1.1% for LRW and 1.3% for Fig. 3. Mean of macrostructural measures derived from T1-w images across all sites across C1-C3 segments of the cervical cord. The mean per subject was calculated across the scan and rescan measurements as well as across sites. Box plot = interquartile range (q3 -q1); bold horizontal line = median, feathers = data range of values across subjects. Cross = outlier, which exceeds the range of (q1-1.5 × [q3 − q1]; q3 + 1.5 × [q3 − q1]), where q1 and q3 are the 25th and 75th percentiles of the sample data, respectively. CSA: spinal cord cross-sectional area, LRW: spinal cord left-right width, APW: spinal cord anterior-posterior width. Fig. 4. Intra-site CoV of parameter maps for individual C1, C2 and C3 segments scanned across all sites. The intra-site CoV was calculated for every subject within each site. Box plot = interquartile range (q3 -q1); bold horizontal line = median, feathers = data range of values across subjects. Cross = outlier, which exceeds the range of (q1-1.5 × [q3 − q1]; q3 + 1.5 × [q3 − q1]), where q1 and q3 are the 25th and 75th percentiles of the sample data, respectively.

Fig. 5.
Inter-site CoV of parameter maps for the C1, C2 and C3 levels of the cervical cord acquired across all sites. The inter-site CoV was calculated for every subject between sites. Box plot = interquartile range (q3 -q1); bold horizontal line = median, feathers = data range of values across subjects. Cross = outlier, which exceeds the range of (q1-1.5 × [q3 − q1]; q3 + 1.5 × [q3 − q1]), where q1 and q3 are the 25th and 75th percentiles of the sample data, respectively.

Relative bias measures
First, we assessed the relative bias of each site separately, which is a measure of the difference between the measurements of one site and the mean of all six sites and thereby allows to evaluate the systemic bias caused by different MRI acquisition protocols at different sites. The RMS of the relative bias was 9.4% for MT, 5.4% for PD, 7.2% for R1, and 12.0% for R2 * ( Fig. 7 ).
Next, we divided the six sites into two separate subgroups with respect to the MRI vendor and compared these subgroups to the whole dataset including all six scanning sites. The sub-groups were a) all sites equipped with siemens scanner and b) all sites equipped with Philips scanner. Of note, the CoVs of most of the MPM maps at the cervical cord were in a similar range comparing the data acquired across the two vendors (Siemens and Philips), except for R2 * maps ( Fig. 8 ).

Discussion
This study determined the intra (scan-rescan) repeatability and intersite comparability of cervical cord MPM data across six different clinical sites equipped with MRI scanners from two different vendors (Siemens and Phillips) by a traveling heads study. The intra-site (scan-rescan) repeatability CoV was between 2.5% and 7.5% across five subjects and inter-site comparability CoV was between 6.1% and 12% across six sites. The intra-and inter site CoV of the cervical cord morphometric measures ranged between 1.1% and 4.0% for CSA, LRW and for APW. By including diverse hardware and software, this study captures the situation across multicenter clinical trials and provides a realistic reference for the performance of MPM protocol in multicenter studies using widely available product sequences and processing methods.

MPM values in the cervical cord
The mean values of MT, PD, R1, and R2 * in cervical levels (C1-C3) across participants were in line with previous reports at cervical cord levels. For instance, the mean value of the R1 and R2 * (mean ± SD: 0.88 ± 0.09 s − 1 and 19.98 ± 3.00 s − 1 respectively) were in the range of reported values in the spinal cord (R1; 87.7 ± 9.04 s − 1 & R2 * ; 20.0 ± 3.0 s − 1 . Likewise, mean PD in the spinal cord was in the range of values reported previously in the cervical spinal cord ( Lévy et al., 2018 ). Of note, the mean MT in this study was 3.28 ± 0.45 p.u., and hence was in line with those previously reported ( Lévy et al.,  Fig. 6. Intra and inter-site CoV of the spinal cord macrostructural measures across C1, C2 and C3 level of cervical cord. Cross-sectional area (CSA), left-right width (LRW) and anterior-posterior width (APW) at cervical level C1-C3 level. Box plot = interquartile range (q3 -q1); bold horizontal line = median, feathers = data range of values across subjects. Cross = outlier, which exceeds the range of (q1-1.5 × [q3 − q1]; q3 + 1.5 × [q3 − q1]), where q1 and q3 are the 25th and 75th percentiles of the sample data, respectively. 2018 ; . However, there is no gold-standard protocol and value for the MT saturation due to the fact that the MT pulse is limited to the vendors' default product sequence settings. Therefore, the MT values were harmonized across vendor platforms in the postprocessing in this study ( Leutritz et al., 2020 ). In general, residual deviations may exist in multi-center studies due to methodological and instrumental differences as well as inter-individual biological variation. For instance, the different approaches in correcting B1 + / − field and incomplete spoiling biases ( Preibisch and Deichmann, 2009 ) may result into slightly different outcomes. Moreover, gradient non-linearity may contribute to some of the differences. Of note, in this study, the different field strengths between the two vendors (Siemens = 2.89T, Philips = 3.00T) were connived which may lead to a small bias of approximately 1.4% in T1 relaxation ( Rooney et al., 2007 ).

Intra-and inter-site reproducibility and relative bias of MPM maps
For the cervical cord, we chose a ROI-based analysis and calculated the reproducibility of averaged measures per vertebral level in the spinal cord which is a common approach in the spinal cord MRI processing ( Samson et al., 2013 ;Martin et al., 2016 ;Lévy et al., 2018 ;Lukas et al., 2021 ). Using this approach, the MPM maps in the cervical cord showed a high intra-site reproducibility in the range of 2.6% and 7.5% across cervical levels C1-C3. Similarly, the inter-site CoV was in a range of 6.5% and 12% for all the MPM maps. Although, our cervical cord data showed higher intra-site repeatability for some maps (R2 * ) compared to the brain data of the same subject presented previously ( Leutritz et al., 2020 ), a direct comparison between CoVs of spinal cord and the brain data is not recommended as Leutritz et al. calculated the reproducibility of brain MPM by applying a voxel-wise approach.
The observed CoV of the spinal cord MPM parameters in this study was 3 times lower (in the range of 2.5% and 12%, Figs. 3 and 4 ) than the trauma-related effect sizes shown in longitudinal studies of spinal cord injury which is in the range of 30% changes using ROI based approach ( Freund et al., 2013 ). Hence, the presented MPM protocol is expected to reliably detect and characterize the trauma-related effects in longitudinal multicenter studies ( Leutritz et al., 2020 ). The root-mean-square (RMS) value (x RMS ) of the inter-site biases were between 5.4% and 12% for different sites for MT, PD and R1 which is comparable in the brain. However, the R2 * variability was still higher than the other MPM parameters both in the spinal cord and the brain ( Leutritz et al., 2020 ), partly due to higher intra-site CoV of R2 * maps and systematic differences between the two MRI manufactures included in this study. The source of higher variability of R2 * could also arise from magnetic field inhomogeneities which may be corrected by the approach suggest by ( Baudrexel et al., 2009 ).

Morphometric measures of the cervical cord
The mean cervical cord cross-sectional area (62.04 ± 4.85 mm 2 ) and the APW and LRW (7.28 ± 0.34 mm, 10.79 ± 0.84 mm, respectively) was calculated based on the T1-weighted scan derived from the MPM protocol and agreed with previous reported values ( Lundell et al., 2011 ;Cohen-Adad et al., 2021b ;Lukas et al., 2021 ). It is known that spinal cord cross-sectional area estimates can systematically differ between the various segmentation software methods ( Prados et al., 2017 ;Weeda et al., 2019 ;Lukas et al., 2021 ). Hence, the values calculated in this study compared to studies using the JIM software ( Horsfield et al., 2010 ) which has implemented a semi-automated algorithm to calculate the CSA is generally bigger (approximately mean 76·7 ± 4·6  mm 2 ) ( Freund et al., 2013 ;Grabher et al., 2017 ;Seif et al., 2019 ;Vallotton et al., 2021 ). Moreover, differently weighted MR images also impact the calculation of the cord area due to different contrasts between the spinal cord and CSF. This speaks for the need of future studies that investigate the effects of different MR contrast on the size of the calculated spinal cord cross-sectional area.

Considerations and limitations
The small sample size of five participants may not represent the general population and its variability. Thus, care should be taken when extrapolating these results to different patient or healthy control groups. However, we believe that most of the study characteristics are fundamental and will be only modulated by the population studied. Furthermore, the results of this study which is based on ROI analysis must be interpreted carefully when comparing with reproducibility studies that applied a different approach such as voxel-based measures, as the reproducibility measures may differ considerably between these approaches given ROI based analysis results into smaller CoVs. Moreover, different MRI scanners are equipped with different MT pulses which resulted in about 20% difference in MT saturation values in this study ( Leutritz et al., 2020 ). Therefore, the MT maps were harmonized during the post-processing by a linear rescaling of the MT values for both brain and spinal cord ( Leutritz et al., 2020 ). This may impact the intersite bias due to different pulse designs. The MT harmonization would only marginally influence the maps. However, rescaling MT values may result in different values compared with previous studies, for example, due to differences between custom-made sequences with optimized MT pulse or different standardized MT sequences.

Conclusion
This study corroborates a high intra-and inter-site reproducibility of the MPM data of the cervical cord obtained from the MPM protocol that covers the whole brain and cervical cord at 3T MRI using Philips and Siemens vendor sequences. The measurements were comparable across sites and scans, as reflected by a low inter-site bias and highly reproducible. Our MPM protocol can be readily applied in qMRI single-and multi-site studies as we only used optimized vendor product sequences for the data acquisition and the open source SCT-toolbox ( De Leener et al., 2017 ) and hMRI toolbox ( Tabelow et al., 2019 ) for the data processing. The results enable and crucially improve the comparability of multicenter studies towards standardized neuroimaging biomarkers applied in clinical trials.

Data and code availability statement
The anonymized data that support the findings of this study are available from the corresponding author upon request.

Study funding
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the grant agree-  (IRP-2022-01-158).

Declaration of Competing Interest
The Max Planck Institute for Human Cognitive and Brain Science has institutional research agreements with Siemens Healthcare. NW was a speaker at an event organized by Siemens Healthcare and was reimbursed for the associated travel expenses.

Data availability
Data will be made available on request.