Simultaneous voxel‐wise analysis of brain and spinal cord morphometry and microstructure within the SPM framework

Abstract To validate a simultaneous analysis tool for the brain and cervical cord embedded in the statistical parametric mapping (SPM) framework, we compared trauma‐induced macro‐ and microstructural changes in spinal cord injury (SCI) patients to controls. The findings were compared with results obtained from existing processing tools that assess the brain and spinal cord separately. A probabilistic brain‐spinal cord template (BSC) was generated using a generative semi‐supervised modelling approach. The template was incorporated into the pre‐processing pipeline of voxel‐based morphometry and voxel‐based quantification analyses in SPM. This approach was validated on T1‐weighted scans and multiparameter maps, by assessing trauma‐induced changes in SCI patients relative to controls and comparing the findings with the outcome from existing analytical tools. Consistency of the MRI measures was assessed using intraclass correlation coefficients (ICC). The SPM approach using the BSC template revealed trauma‐induced changes across the sensorimotor system in the cord and brain in SCI patients. These changes were confirmed with established approaches covering brain or cord, separately. The ICC in the brain was high within regions of interest, such as the sensorimotor cortices, corticospinal tracts and thalamus. The simultaneous voxel‐wise analysis of brain and cervical spinal cord was performed in a unique SPM‐based framework incorporating pre‐processing and statistical analysis in the same environment. Validation based on a SCI cohort demonstrated that the new processing approach based on the brain and cord is comparable to available processing tools, while offering the advantage of performing the analysis simultaneously across the neuraxis.

However, there is no solution that enables the simultaneous assessment of changes occurring in the brain and cord. As a result, most imaging studies fail to implement and analyse the interactions between such remote areas across the central nervous system. Providing a tool that could assess the sequelae of a focal central nervous system (CNS) injury across the entire neuroaxis simultaneously would hold great potential to better understand the temporally and spatially distributed pathophysiological changes (Freund et al., 2016).
Anatomical templates are fundamental to group level comparisons because they naturally define a common space in which analyses can be performed. Templates typically either take the form of mean intensity images or of tissue probability maps (TPMs); the latter can further serve as anatomical priors in segmentation algorithms. A wide variety of brain templates have been developed, serving different purposes, but spinal cord templates are less common. Recently, a template that covers brainstem and spinal cord has been proposed (De Leener et al., 2018), using spinal cord specific tools that cannot be applied to the brain, so the conjoint analyses of the brain and cord is not possible. A principled way of constructing a template is to consider it as an unknown variable in a generative model of large MR datasets (Blaiotta, Freund, Cardoso, & Ashburner, 2018). In this context, a brain-neck template has naturally been obtained by fitting the generative model to collections of images of the brain and spinal cord (Blaiotta et al., 2018).
The present study implements a brain-neck template (covering brain and cervical spinal cord) and validates its use in the context of VBM (Ashburner & Friston, 2005) and VBQ  analyses as implemented in the SPM framework. The extended SPM framework for brain and spinal cord analyses, named 'SPM-BSC', offers the advantage of performing pre-processing and statistics simultaneously in the brain and cervical spinal cord. In addition to voxel-wise analyses in standard (group) space, cord-specific metrics were also computed in native space. These include spinal cord area (SCA), anterior-posterior width (APW) and left-right width (LRW).
Furthermore, mean values across the whole cervical spinal cord of the magnetization transfer saturation (MT), longitudinal relaxation rate (R1 = 1/T1) and effective transverse relaxation rate (R2* = 1/T2*) were also computed. The SPM-BSC analysis was validated by assessing volume and microstructural changes in the brain and spinal cord in SCI-patients relative to healthy controls and comparing the results with pipelines separately optimised for the brain (SPM-BO for 'brain only'), or from the SCT covering only the spinal cord.  (Kirshblum et al., 2011) for motor, light touch, and pinprick score; the Spinal Cord Independence Measure (SCIM) (Catz et al., 2007); and additionally for tetraplegic patients the Graded Redefined Assessment of Strength, Sensibility (Table 1).

| Image acquisition
Participants underwent a 3D T1-weighted MPRAGE (magnetizationprepared rapid acquisition gradient echo) scan (whole-brain including the cervical cord down to C5 level) on a 3 Tesla MRI scanner (Magnetom Verio [n = 35 subjects] or SkyraFit [n = 18 subjects] Siemens Healthcare, Erlangen, Germany) with the following parameters: field of view (FOV) of 224 × 256 × 176 mm 3 , isotropic resolution of 1 mm 3 , repetition time (TR) 2,420 ms, Echo time (TE) 4.18 ms, flip angle (α) 9 , inversion time 960 ms, and readout bandwidth of 150 Hz per pixel. The system was equipped with a 16-channel radiofrequency (RF) receive head and neck coil and RF body transmit coil.
To assess microstructural changes, a multi-parameter mapping (MPM) protocol, based on multi-echo 3D FLASH sequences (Leutritz et al., 2020;Tabelow et al., 2019), was performed. These scans covered the whole brain and cervical spinal cord down to level C4 with 1 mm isotropic resolution and a field of view of 240 × 256 mm 2 (matrix size of 240 × 256) with 176 partitions in a total scan time of 23 min. Parallel imaging with a speed up factor of 2 was used in the phase-encoding direction (anterior-posterior) using a generalised auto-calibration partially parallel acquisition algorithm (GRAPPA) and a partial Fourier acquisition with a 6/8 sampling factor was used in the partition direction (left-right) with a readout bandwidth of 480 Hz/pixel. Each set of echoes was acquired using a different TR and flip angle (α) to achieve images with either T1-weighting: 25 ms/23 , PD-weighting: 25 ms/4 , or MT-weighting: 37 ms/9 with off-resonance RF pulse (Gaussian shape with pulse length of 10 ms, flip angle of 500 , frequency 1.2 kHz and bandwidth of 192 Hz) prior to excitation. Echoes were acquired at six equidistant echo times (TE) from 2.46 ms to 17.22 ms for all weightings, with an additional echo at 19.68 ms for the PD-weighted and T1-weighted volumes.

| Image analysis
VBM and VBQ analyses in SPM12 were extended to include cervical spinal cord ('SPM-BSC') by incorporating a new probabilistic atlas of the brain and neck (Blaiotta et al., 2018). Validation of this approach was performed by analysing a cohort of 53 individuals comprised of 30 patients with chronic traumatic SCI and 23 healthy participants. The SPM-BSC analysis was used to assess volume and microstructural differences within the brains and spinal cords of the SCI-patients relative to the healthy controls. Results obtained with the proposed SPM-BSC analysis were compared with those obtained using existing tools that separately assess either the brain or spinal cord only. Therefore, MRI data were analysed using: • SPM-BSC: the proposed VBM/VBQ pipeline implemented in SPM12 using a brain and spinal cord template enabling statistical analysis across the entire neuroaxis.
• SPM-BO: the established VBM/VBQ pipelines available in SPM12 within the brain only.

| Simultaneous analysis of the brain and spinal cord (SPM-BSC)
The template incorporating brain and cervical spinal cord (named brain-neck template) had been computed (Blaiotta et al., 2018) using a generative semi-supervised modelling approach (Blaiotta, Jorge Cardoso, & Ashburner, 2016). The algorithm is based on a Gaussian mixture model of MR intensities, where the intensity distribution within each tissue type is assumed Gaussian, and the template encodes prior probabilities of belonging to each tissue (Blaiotta et al., 2016). A total of 12 weakly supervised components were optimised and, after algorithm convergence, they were combined into seven tissue maps based on visual inspection. In particular, a single component was used to model each of grey matter (GM), white matter (WM), cerebrospinal fluid (CSF) and fat; two components were used to model non-neural tissues; three to model soft tissues and three to model a mixture of bone and air. Finally, the template was registered to the existing SPM12 tissue probability template using a 12-parameter affine registration. This template covers the brain and spinal cord down to level C3.

| SPM-BSC analysis in MNI space
The brain-neck template (Blaiotta et al., 2018) was used to segment the T1-weighted structural MPRAGE images of each individual using the unified segmentation algorithm (Ashburner & Friston, 2005) in SPM12. The algorithm can use more than one Gaussian distribution to model intensities within each tissue class. Here, one Gaussian was specified for GM, WM, CSF and fat; two Gaussians were specified for non-neural tissues and three Gaussians were specified for each of the other two tissue classes (soft tissues and a mixture of bone and air).
To assess morphological changes in the spinal cord, the native-space GM and WM tissue maps were combined to form a neural tissue (NT) class. Then, the GM, WM and NT maps were spatially normalised to MNI space with Dartel (Ashburner, 2007), and modulated by the Jacobian determinants of the deformations (Good et al., 2001). Finally, an isotropic Gaussian kernel of 6 mm full width at half maximum (FWHM) was applied to the modulated tissue maps ( Figure 1a). The total intracranial volume (TIV) was computed from the sum of the GM, WM, and CSF volumes (Ridgway, Barnes, Pepple, & Fox, 2011).
The VBQ analysis (Draganski et al., 2011;Tabelow et al., 2019) was based on the quantitative metrics of MT, R1 and R2* derived from the MPM using the in vivo histology MRI (hMRI) toolbox (Tabelow et al., 2019). Inhomogeneity of the RF transmit field was corrected using UNICORT . The VBQ analysis included segmentation of the MT map into GM, WM and CSF, covering brain and spinal cord, using the unified segmentation approach (Ashburner & Friston, 2005) and the brain-neck template (as used when segmenting the MPRAGE for the VBM analysis). Segmenting the MT map ensures that the MPM maps and the tissue probabilities are aligned in native space. To assess microstructural changes in the spinal cord, NT maps were computed by again combining the GM and WM tissue maps from MT segmentation in native space. All tissue maps were transformed to MNI space using the Dartel algorithm (Ashburner, 2007). MPMs were warped to the MNI space using participant-specific warping fields generated from the Dartel algorithm and finally smoothed using an isotropic Gaussian kernel filter with 6 mm FWHM using the VBQ smoothing approach with the corresponding tissue probability map (Draganski et al., 2011) ( Figure 1b).
Smoothing the spinal cord MRI leads to undesired partial volume effects with CSF. Therefore, we investigated different smoothing kernels to identify the best scale space (Worsley, Marrett, Neelin, & Evans, 1996), for spinal cord analysis. We applied isotropic smoothing kernels with FWHM between zero (i.e., no smoothing) and 6 mm, with step of 1 mm, prior to the statistical analysis.

| SPM-BSC analysis in subject space
To compute morphometric metrics in the spinal cord, such as SCA, APW and LRW, the native-space NT tissue probability maps were thresholded at 0.5 to generate a binary mask of spinal cord. Next, an ellipse was fitted to the boundary of the mask in each transverse slice and the corresponding area (SCA) and main cord axes (LRW and APW) were extracted. Microstructural metrics were computed as the mean of the MT, R1 and R2* values within the masked region of the NT probability maps (computed from the MT map segmentation) covering cervical segments C1 to C3.

| Brain analysis (SPM-BO) in MNI space
MPRAGE anatomical images were analysed using the standard VBM approach, which segments and aligns the brain but does not account for the spinal cord. This procedure involved segmentation into GM, WM and CSF using unified segmentation (Ashburner & Friston, 2005) with the tissue priors released with the SPM12 software. The GM and WM tissue probability maps were subsequently warped into standard MNI space using Dartel, before being smoothed with an isotropic Gaussian kernel with 6 mm FWHM. In order to provide voxel-wise statistic in PAM50 space using SCT, the smoothed MPM maps were warped into this template space using the warping fields derived from the segmentation of the MT maps.

| Regions of interest
An ROI approach was applied based on the literature findings (Freund et al., 2013;Grabher et al., 2015). This was necessary to be equally sensitive to brain and cord changes using data from the different methods. The following ROIs were used in the voxel-wise analysis of the brain: bilateral motor cortex M1 (precentral gyrus), bilateral somatosensory S1 cortices (postcentral gyrus), thalamus and corticospinal tract, which were extracted using the SPM Anatomy toolbox (Eickhoff et al., 2007 According to the ICC guidelines (Koo & Li, 2016), ICC values can be partitioned to indicate poor (less than 0.5), moderate (between 0.5 and 0.75), good (between 0.75 and 0.9), and excellent (greater than 0.90) reliability between methods.

| Associations between spinal cord MPM maps and clinical outcomes
Spinal cord MRI data in the template spaces (from both SPM-BSC and SCT analyses separately) were used to assess the associations between MPM measures and clinical outcomes (LEMS, light-touch, pinprick and SCIM scores), using regression models in SPM12. All tests were corrected for age, gender, scanner and total intracranial volume (TIV) as covariates of no interest and family-wise error (FWE) correction was applied to account for multiple comparisons (Friston et al., 1994) using a threshold of p = 0.05 at the peak-level.  Table 2.

| Structural changes in the brain using SPM-BO
The ROI-based analysis revealed statistically significant atrophy (zscore = 4.2, x = 2, y = −11, z = 5, p = .011) and myelin-sensitive R1 reduction (  F I G U R E 2 Overlay of statistical parametric maps (uncorrected p < .001, for illustrative purposes) showing morphometric and microstructural changes in SCI-patients compared to healthy controls using data pre-processed via the simultaneous brain and spinal cord (SPM-BSC) approach in (a) and the standard brain approach in (b). Both pre-processing methods showed myelin-sensitive R1 reduction in the thalamus (yellow); myelinsensitive R1 and MT reduction in the bilateral sensory-motor cortex (R1, yellow; MT, red); and myelin-sensitive R1 reduction in right corticospinal tract in SCI-patients compared to healthy controls (R1, yellow; MT, red). In addition, the SPM-BSC approach (a) showed atrophy (green) and myelin-sensitive MT reduction in the cervical spinal cord of SCI-patients compared to healthy controls. The colour bar indicates the t score

| Spinal cord changes using voxel-wise analysis in template space
Data pre-processed using the SPM-BSC approach showed spinal cord atrophy (from the VBM analysis) and decreased MT (form the VBQ analysis) in SCI-patients compared to controls for all smoothing sizes (Table 4). Decreased spinal cord R1 was found in SCI-patients compared to healthy controls only if no smoothing was applied (z-score = 3.82, x = 0, y = −49.5, z = −78, p = .042, in brain-neck template space) (Figure 3-a2).

| Spinal cord changes using extracted metrics in subject space
Morphometric analysis from data pre-processed with the SPM-BSC approach, revealed reduced SCA and APW in SCI-patients compared to healthy controls (p < .05); while the LRW showed only a trend (p = .08) in the cervical spinal cord. Similarly, data pre-processed with the SCT approach revealed reduced SCA and APW (p < .05) in SCIpatients compared to healthy controls; while the LRW showed only a trend (p = .06) (Figure 4a and Table 5). Microstructural analysis from data pre-processed with the SPM-BSC approach, revealed reduced myelin-sensitive MT and R1 (p < .05) in SCI-patients compared to healthy controls; whereas measures of R2* in the cervical cord were not found to be different (p = .49). Similarly, data pre-processed with T A B L E 2 Results from region of interest (ROI) analysis using voxel-based morphometry (VBM) and voxel based quantification (VBQ) from 'Brain-neck SPM' and 'Brain only' methods ROI z-score the SCT approach revealed reduced myelin-sensitive MT (p < .001) in SCI-patients compared to healthy controls, but no differences were found in either the R2* maps (p = .55) or the R1 maps (p = .62) ( Figure 4b and Table 5).

| DISCUSSION
In this study, a template covering the brain and the cervical spinal cord was validated within the SPM12 framework. This was done by comparing the results from the between-group effect of SCI patients and healthy controls on qMRI data in the brain and in the spinal cord using the new SPM-BSC analysis and established SPM brain and spinal cord toolbox analysis. In addition, good to excellent ICCs were observed in the majority of brain ROIs. Thus, the proposed SPM-BSC approach enables analysis of equivalent data, sensitive to disease-associated effects as brain or spinal cord, as established region-specific analysis tools but with the advantage of being embedded within the simple unified SPM framework.
SCI patients exhibited lower R1 values in the thalamus, corticospinal tract and sensory-motor cortex compared to healthy controls. Patients also displayed lower MT values than healthy controls in the sensory-motor cortex. These changes are consistently found across analysis methods (see Table 2) and are in line with previous studies (Grabher et al., 2015;Seif et al., 2018), supporting the conclusion that SPM-BSC approach is as sensitive as SPM-BO for detecting group differences.
In addition, the ICC showed good to excellent values for most of the brain ROIs. Moderate ICCs were reported in the thalamus and corticospinal tract for MT and volume maps. These lower values might be due to discrepancies between the two templates, which were visually assessed by examining the difference between the standard and brain-neck templates. The observed discrepancies can be attributed to the alignment strategy; in particular the standard brain template is smoother than the brain-neck template since the latter used more precise alignment during the template generation. More direct comparison would be aided by using the same alignment strategy during the generation of the template and subsequent segmentation for both analysis approaches.
Voxel-wise analysis in the spinal cord using the SPM-BSC approach showed that SCI patients exhibited lower MT values than healthy controls, regardless of the size of the smoothing kernel used.
However, R1 maps did not show the same results over all smoothing sizes and reduced R1 values were found only for unsmoothed data ( Figure 3). This suggests that only focal differences were present in the R1 map within the cord and were obscured by partial volume effects resulting from the smoothing process. The SCT analysis using the PAM50 template showed reductions in MT and R1 values in SCIpatients compared to healthy controls; and similar findings were F I G U R E 3 Overlay of statistical parametric maps (uncorrected p < .001, for illustrative purposes) showing statistically significant differences in SCI-patients compared to healthy controls in the cervical spinal cord using different smoothing sizes and data analysed with SPM-BSC (in a) and SCT (in b) approaches. The colour bar indicates the t score found using smoothed data. The discrepancy between findings of reduced R1 between the two approaches in the smoothed case might be due to the smoothing algorithm applied: in the SPM-BSC analysis, weighted isotropic Gaussian kernel smoothing was applied within the NT tissue class using the same VBQ approach applied in the brain (Draganski et al., 2011), enabling simultaneous assessment of the brain and neck using random field theory (Eklund, Nichols, & Knutsson, 2016). In contrast, in the SCT approach, a centerline smoothing algorithm was used, which works by straightening the cord, applying a 1D Gaussian smoothing kernel along its length and then un-straightening the cord back into the original space (De Leener et al., 2017). The 1D nature of this latter approach may be more robust to partial volume effects, even within the cord, than the isotropic smoothing adopted in the VBQ approach. Nonetheless, both tools showed myelin sensitive MT changes that were consistent for all smoothing sizes. Further investigation of the impact of smoothing on voxel-wise analyses in the spinal cord is needed.
Spinal cord native-space analysis revealed atrophy and microstructural changes in SCI-patients compared to healthy controls. In particular, both approaches (SPM-BSC and SCT) showed atrophy in the spinal cord, more specifically in APW, and myelin-sensitive MT reduction in SCI-patients compared to healthy controls. However, different results were obtained for R1, where the SPM-BSC approach showed differences between SCI-patients compared to healthy controls, whereas in SCT this difference was not statistically significant.
However, results from the analysis in PAM50 space (SCT) were more in line with those from the SPM-BSC approach. Atrophy and microstructural changes in the spinal cord observed with the SPM-BSC and SCT approaches in patients with SCI are in line those in the literature (Freund et al., 2019). This demonstrates that the SPM-BSC approach can detect the same injury-related differences between SCI-patients and healthy controls as existing tools that are dedicated to either brain or cord. This is achieved without a loss of sensitivity and is the case for both native-space spinal cord metrics and template-space voxel-wise analyses.
The association analysis between spinal cord microstructure and clinical disability showed that higher MT values in the spinal cord were associated with better sensory and motor outcome. The F I G U R E 4 Structural changes in the cervical spinal cord applying SPM-BSC and SCT tools: (a) morphometric analysis showing reduced crosssectional area and anterior posterior width in SCI-patients compared to healthy controls (p < .05). (b) Microstructural analysis showing reduced MT in SCI-patients compared to healthy controls (p < .05). The * indicates statistically significant differences (p < .05) between the connected groups SPM-BSC approach was able to replicate the associations between spinal cord microstructure and functional impairment in SCI-patients that were identified by the SCT analysis, though the latter did show stronger correlation (higher R values). This might again be due to the different smoothing approaches used by the two methods, that is, SCT uses an adaptive Gaussian kernel oriented along the spinal cord centreline (De Leener et al., 2017). There are no standards for the correct shape and size of smoothing kernel to be used in the spinal cord, warranting further investigation on this topic.
Several limitations must be considered when using the new SPM-BSC approach. First, the brain-neck template covers only cervical spinal cord until level C3. However, a template spanning more spinal cord levels could be generated using the same algorithm (Blaiotta et al., 2018) by including MRI data covering more spinal cord segments. This will be the focus of future work. Furthermore, the GM/WM classification in the spinal cord is sub-optimal. As a result, the probability maps of GM and WM were combined in the present study to generate a probability map of NT, which is a nonspecific measure of spinal cord volume. In future, better GM/WM classification within the spinal cord might be achieved by combining ultra-high resolution spinal cord data and high-resolution brain MRI during template creation. Lastly, results from data preprocessed with the standard brain-based SPM approach showed slightly higher z-scores (see Table 3), which could be due to differences in alignment due to the templates and segmentation inputs used. F I G U R E 5 Associations between clinical measures and spinal cord MT maps pre-processed using SPM-BSC (in a) and SCT (in b); input data were smoothed with a 3 mm isotropic Gaussian kernel for SPM-BSC and with sigma 3 mm for SCT approach using SCT tools