Relationship between brainstem neurodegeneration and clinical impairment in traumatic spinal cord injury

Background Brainstem networks are pivotal in sensory and motor function and in recovery following experimental spinal cord injury (SCI). Objective To quantify neurodegeneration and its relation to clinical impairment in major brainstem pathways and nuclei in traumatic SCI. Methods Quantitative MRI data of 30 chronic traumatic SCI patients (15 with tetraplegia and 15 with paraplegia) and 23 controls were acquired. Patients underwent a full neurological examination. We calculated quantitative myelin-sensitive (magnetisation transfer saturation (MT) and longitudinal relaxation rate (R1)) and iron-sensitive (effective transverse relaxation rate (R2*)) maps. We constructed brainstem tissue templates using a multivariate Gaussian mixture model and assessed volume loss, myelin reductions, and iron accumulation across the brainstem pathways (e.g. corticospinal tracts (CSTs) and medial lemniscus), and nuclei (e.g. red nucleus and periaqueductal grey (PAG)). The relationship between structural changes and clinical impairment were assessed using regression analysis. Results Volume loss was detected in the CSTs and in the medial lemniscus. Myelin-sensitive MT and R1 were reduced in the PAG, the CSTs, the dorsal medulla and pons. No iron-sensitive changes in R2* were detected. Lower pinprick score related to more myelin reductions in the PAG, whereas lower functional independence was related to more myelin reductions in the vestibular and pontine nuclei. Conclusion Neurodegeneration, indicated by volume loss and myelin reductions, is evident in major brainstem pathways and nuclei following traumatic SCI; the magnitude of these changes relating to clinical impairment. Thus, quantitative MRI protocols offer new targets, which may be used as neuroimaging biomarkers in treatment trials.


Introduction
Traumatic spinal cord injury (SCI) is a devastating condition and causes permanent sensorimotor loss and autonomic dysfunction in most patients, with no cure currently available (Dietz and Fouad, 2014). Usually patients show some degree of recovery which levels off within two years after injury. Using computational neuroimaging approaches, rapid and dynamic trajectories of neurodegenerative processes have been identified above the level of injury that accompanied the recovery. Crucially, the magnitude of neurodegeneration was associated with clinical impairment (Freund et al., 2013;Grabher et al., 2015).
Besides neurodegeneration at the spinal and cortical level (Beaud et al., 2008;Jirjis et al., 2015), retrograde and transneuronal degeneration has been shown in experimental SCI in brainstem pathways (Jirjis et al., 2015;Jones and Pons, 1998) and nuclei (Jones and Pons, 1998;Kwon et al., 2002;Wannier-Morino et al., 2008). The brainstem is phylogenetically highly conserved in mammals and plays a key role in motor (Lemon, 2008) and sensory function (Benarroch, 2012;Liao et al., 2015). Important substructures of the motor system entail the rubrospinal system (i.e. execution of precise limb movements), the vestibulospinal system (i.e. balance and posture), the reticular formation (i.e. initiates and coordinates limb movements and postural support), and the corticospinal system (i.e. skilled motor function) (Lemon, 2008), while the dorsal column nuclei and medial lemniscus (Liao et al., 2015) and the periaqueductal grey (PAG) (Benarroch, 2012) are involved in sensory processing and pain modulation. Crucially, structural reorganization of brainstem pathways and nuclei has been associated with functional recovery following experimental SCI (Zaaimi et al., 2012;Zörner et al., 2014). Thus, understanding trauma-induced pathophysiological processes affecting the brainstem pathways and nuclei might offer crucial insights into neurodegeneration and plasticity.
However, the brainstem is still understudied in human SCI as accurate and sensitive neuroimaging tools targeting the brainstem have only recently become available (Lambert et al., 2013b). First attempts using neuroimaging approaches provided evidence of brainstem atrophy (i.e. volume loss) (Freund et al., 2013(Freund et al., , 2012Grabher et al., 2015;Wrigley et al., 2009) and plasticity (i.e. volume increases) during intensive training (Villiger et al., 2015) in human SCI. Recent improvements in quantitative MRI (qMRI) techniques now allow quantification of the underlying microstructural changes  and segmentation of individual brainstem pathways and nuclei (Lambert et al., 2013b). This is possible because different MR contrasts can be used to calculate quantitative maps (magnetisation transfer saturation (MT), longitudinal relaxation rate (R1), effective transverse relaxation rate (R2*)), which are sensitive to myelin (Schmierer et al., 2004;Turati et al., 2015) and iron (Stüber et al., 2014). Such maps can be used for multiparametric brainstem tissue segmentation (Lambert et al., 2013b). Myelin reductions have been shown to accompany atrophic changes in the cord and cortex, thus offering complementary insights into the sequela of SCI (Freund et al., 2013;Grabher et al., 2015). Furthermore, iron accumulation due to myelin breakdown has been reported in SCI (Kroner et al., 2014;Sauerbeck et al., 2013).
Here, we combined voxel-based quantification and multiparametric tissue segmentation to address our hypotheses that after traumatic chronic SCI, (1) atrophy and myelin reduction are evident in major brainstem pathways and nuclei and, (2) that the extent of atrophy, myelin reduction and iron accumulation relates to clinical impairment, lesion level and severity.

Participants and study design
We recruited 30 individuals with a chronic traumatic SCI (3 female) and 23 healthy participants (10 female) at the University Hospital Balgrist between August 2011 and May 2015. Fifteen patients were tetraplegic and fifteen paraplegic. All patients were treated surgically for decompression. No participant reported a history of medical, neurological, or psychiatric disorders and all were eligible for MRI examinations.
To define the level of sensory and motor impairment, the most caudally intact dermatome for light touch and pinprick sensation (2/2 points) and motor function were considered, respectively (according to the ISNCSCI protocol). Lesion-level (neurological level of injury) was defined as the most caudal segment of the cord with intact sensation and motor function against gravity (min. 3/5 points), provided that motor and sensory function above this segment were normal. Lesion completeness was defined as having no motor and sensory function preserved in the sacral levels S4/5 (AIS A).
All participants gave informed written consent prior to study enrolment. The study protocol was in accordance with the Declaration of Helsinki and approved by the Ethics Committee of the Canton Zurich (reference number: EK-2010-0271).

Image acquisition
All participants' structural whole-brain data, including the cervical cord up to vertebra C5, were acquired on a 3T Magnetom MRI scanner (Siemens Healthcare, Erlangen, Germany). The system was equipped with a 16-channel radiofrequency (RF) receive head and neck coil and RF body transmit coil. A multiecho 3D FLASH (fast low-angle shot) sequence, with the following parameters, was used within a wholebrain multiparameter mapping (MPM) qMRI protocol (Draganski et al., 2011;Weiskopf et al., 2011): field of view (FoV) of 240 × 256 mm 2 , matrix size 240 × 256, isotropic resolution of 1 mm, GRAPPA parallel imaging in phase-encoding direction (anterior-posterior) with speed-up factor of 2, partial Fourier acquisition with 6/8 sampling factor in partition direction (left-right), and a readout bandwith of 480 Hz per pixel. Different weightings were predominantly achieved by choosing repetition time (TR) and flip angle (α): (1) T1-weighted (T1w): 25 ms/ 23°, (2) proton density-weighted (PDw): 25 ms/4°, and (3) MTweighted (MTw): 37 ms/9°with off-resonance RF pulse prior to excitation. Echoes were acquired at seven equidistant echo times (TE) from 2.46 ms to 17.22 ms for all volumes, with an additional echo at 19.68 ms for PDw and T1w.

Image pre-processing
The acquired T1w, PDw, and MTw echoes were first averaged to increase the signal to noise ratio (SNR) and then used to calculate quantitative maps of MT and R1 (Draganski et al., 2011;Weiskopf et al., 2011) in MATLAB (MathWorks, Natick, MA). R2* was calculated from the log signal of the PDw echoes. UNICORT  was used to correct RF transmit field inhomogeneity.

Brainstem template generation
We first generated brainstem tissue probability maps (TPMs) for the spatial alignment of brainstem sub-structures in our study cohort and to increase sensitivity for pathophysiological processes. Before generating the TPMs, we extracted the brainstem from quantitative maps from a longitudinal qMRI dataset of 29 subjects over four time points (Freund et al., 2013;Grabher et al., 2015) by label propagation using a set of brain labels (Neuromorphometrics Inc., Somerville, USA). Subsequently, whole-brain deformation fields were derived by segmenting the MT maps (Ashburner and Friston, 2005) and then applying a diffeomorphic image registration algorithm (Ashburner, 2007). The derived deformation fields enabled the extracted qMRI brainstem data to be transformed to the MNI space.
We then used a multivariate Gaussian mixture model to generate brainstem TPMs (Hasselblad, 1966). Such a model assumes that the observed image intensities are drawn from a set of multivariate Gaussian probability density functions, where each Gaussian captures the intensity distribution of one single tissue type. Additionally, we introduced locally-varying, unknown tissue priors, which are learned directly from the observed data, thus providing a set of populationspecific, average-shaped TPMs (Blaiotta et al., 2016;Lambert et al., 2013b). The statistical Gaussian mixture model was fit to the spatially normalized qMRI brainstem data, using the Expectation-Maximization algorithm (Moon, 1996), which is a general and well-established technique to obtain maximum likelihood or maximum a posteriori estimates of the model parameters, for probabilistic latent variable models. Within the neuroimaging community, such a method has been extensively validated for the classification of neural tissue types from MR data (Ashburner and Friston, 2005;Blaiotta et al., 2016;Lambert et al., 2013b). The resulting seven brainstem TPMs (classes 1-7) are shown in Fig. 1 and contained, amongst others, the red nucleus (RN) (class 6), cerebral crus including the corticospinal tracts (CSTs) (class 6), and PAG (class 3). Anatomical locations were validated using a high-field MRI brainstem atlas (Naidich et al., 2009). The tissue probability maps were subsequently aligned and merged with the whole-brain TPMs provided with SPM12 (http://www.fil.ion.ucl.ac.uk/spm/) (subsequently referred to as modified TPMs), so as to allow a more accurate alignment of the brainstem tissue maps with the individual scans, during the following processing steps. In fact, information derived from the tissues surrounding the brainstem (i.e. grey and white matter) can be effectively used to drive the registration of the individual volumes to the population mean, therefore ensuring more accurate segmentation results.

Voxel-based macrostructural and microstructural analysis of the brainstem
We used our modified TPMs to segment the brain data (using MT and PDw data) of our study population (30 SCI and 23 controls) into grey matter, white matter, cerebrospinal fluid, plus the seven brainstem tissues for each subject (Ashburner and Friston, 2005). Then, a geodesic shooting registration algorithm (Ashburner and Friston, 2011) was used to align the seven brainstem tissues of all subjects and to create a common study population mean (Fig. 2). The estimated deformation fields were used both to compute Jacobian determinant maps for tensor-based morphometry (TBM), and to warp the quantitative maps of MT, R1, and R2* into the study population mean space for voxel-based quantification (VBQ) (Draganski et al., 2011). Due to the lack of gyrification of the brainstem and due to the highly accurate warping algorithm, no smoothing was applied to achieve higher spatial accuracy for small brainstem structures (Lambert et al., 2013a(Lambert et al., , 2013b. Note that using unsmoothed data reduces sensitivity (i.e. fewer false positive results) and increases specificity due to Random Field theory being overly conservative at low smoothness (Nichols, 2012).
We used t-tests within the framework of the general linear model (GLM) to assess morphometric and microstructural differences between individuals with SCI and healthy controls, between tetraplegic and paraplegic patients, and between patients with complete and incomplete lesions. We used regression models to assess the relationship between morphometric and microstructural measures and neurological and functional impairment (AIS, lesion level, LEMS; UEMS, LT, PP, SCIM). Covariates of no interest included age, total intracranial volume and scanner upgrade to control for confounding linear effects in all GLMs (Barnes et al., 2010). Cluster-inference was performed using a cluster-defining threshold of p = 0.001 and a family-wise error (FWE) corrected threshold of p = 0.05 using Gaussian Random Field theory to account for multiple comparisons (Friston et al., 1994) within regions of interest (ROIs) derived from the seven brainstem TPMs. Only significant results (p < 0.05) corrected for FWE are reported. The ROIs (i.e. brainstem TPMs) were used to increase sensitivity for Fig. 1. Brainstem tissue probability maps (classes 1-7). Seven within-brainstem tissue classes were derived from multiparametric brainstem segmentation using a multivariate mixture of Gaussians. They contain brainstem nuclei including the substantia nigra (class 1), the periaqueductal grey (class 3), and the red nucleus and cerebral crus (class 6).

Fig. 2.
Overlay of study population mean onto highresolution Duvernoy Atlas of the brainstem (Naidich et al., 2009). The population mean was derived by registration of the seven brainstem tissues of all subjects and are in correspondence to the obtained brainstem tissue probability maps (see Fig. 1). pathophysiological processes and specificity for anatomical locations in brainstem sub-structures (e.g. CST, RN, PAG).

Volume loss and myelin reductions in brainstem pathways and nuclei
Voxel-wise analysis revealed significant atrophy and myelin reductions in patients compared to healthy controls within the brainstem (Fig. 3, Table 2). Volume loss was observed in the CSTs and medial lemniscus at the level of the medulla (p = 0.017). Lower myelin-sensitive MT was evident in the left CST at the level of the medulla (p = 0.039) and within the PAG (p = 0.001). Lower myelin-sensitive R1 was observed in the CSTs at the level of the medulla (cluster 1: p = 0.003; cluster 2: p = 0.025) and bilaterally in the dorsal medulla (p < 0.001), in the dorsal pontomedullary junction (cluster 1: p = 0.020; cluster 2: p = 0.038) and in the dorsal pons (7 clusters, Table 2). Iron-sensitive R2* did not reveal any significant changes in patients compared to controls. No between group differences were observed in volume and quantitative maps in tetraplegic and paraplegic patients and in patients with complete (i.e. AIS A) and incomplete lesions (i.e. AIS B-E).

Discussion
This study revealed atrophy and myelin reductions (i.e. MT and R1) within major brainstem pathways and nuclei involved in motor and sensory (dys-) function in chronic traumatic SCI. Interestingly, atrophy was observed only in the motor and sensory pathways, whereas myelin reductions occurred also in areas containing brainstem nuclei. Crucially, the magnitude of myelin reductions was related to the extent of motor and sensory impairment. Therefore, these structural alterations in the brainstem could be considered, if reproduced in longitudinal studies, as new targets to monitor impairment and complement assessments in clinical trials following SCI.  retrograde fibre degeneration (Beaud et al., 2008). Interestingly, we did not observe any macrostructural or microstructural changes in the rubrospinal system (i.e. RN). Phylogenetically, the formation of direct cortico-motoneuronal connections (i.e. CST) reaching as far as the lumbar enlargement (Lemon, 2008) has rendered the rubrospinal system less important for the control of movements in man, where fibres terminate in the high cervical cord (Nathan and Smith, 1982). Conversely after experimental SCI, the RN and the rubrospinal fibres, which reach more caudally into the spinal cord, show signs of neurodegeneration (Kwon et al., 2002;Wannier-Morino et al., 2008). Thus, these findings might illustrate the more dominant role of the CST in processes of neurodegeneration and plasticity in the context of functional recovery after human SCI (Dietz and Fouad, 2014). We further observed myelin reductions in areas involved in motor control such as the cerebellar peduncles and ventromedial brainstem pathways (Lemon, 2008;Moulton et al., 2010), independent of atrophy. Reduced postural stability and increased risk of fall in incomplete SCI patients might result from extrapyramidal changes in the vestibular and reticular system (Liechti et al., 2008). Crucially, myelin reductions in the medulla (i.e. vestibular nucleus) and pons (i.e. pontine nucleus) were related to lower functional independence (i.e. SCIM score). Thus we provide evidence that the integrity of the pyramidal and extrapyramidal systemscrucial for postural control and movement coordination and hence functional independenceare disturbed in human SCI. Interestingly, no relationships between clinical impairments and atrophy were observed. Volume loss is rather unspecific to the underlying pathological processes and points to the sensitivity and accuracy of quantitative markers of the myeloarchitecture. Therefore, the extrapyramidal system, in combination with advanced neuroimaging tools, might offer new treatment targets, which could help to increase levels of independence in activities of daily living by improving postural stability (Liechti et al., 2008).

Sensory system neurodegeneration
Next to CST atrophy, we identified atrophic changes in the medial lemniscus. Although we did not observe any concomitant myelin changes, this is suggestive of anterograde fibre degeneration of sensory pathways (Jones and Pons, 1998). The PAG is part of the endogenous pain inhibition system and involved in motor function (Benarroch, 2012). We observed myelin reductions in the PAG which were directly related to impaired pinprick sensation (i.e. pain), a modality that Fig. 3. Volume loss and myelin reductions in brainstem pathways and nuclei in chronic SCI. (A & B) Three dimensional illustration of atrophy (Jacobians, blue) and myelin reductions (R1: yellow, MT: red) in the brainstem. Please note that for illustrative purposes the statistically significant clusters were smoothed with a Gaussian kernel with 1 mm full width at half maximum. Overlay of statistical parametric maps showing atrophy in the corticospinal tracts and medial lemniscus (E) and myelin reductions in the periaqueductal grey (C), in the dorsal pons (D), and in the dorsal medulla (E). conveys afferent information flow via the spinothalamic tract to second order neurons within the PAG. Altered information flow via the nociceptive dorsal horn neurons and the spinothalamic tract (Haefeli et al., 2013) into the PAG might therefore contribute to the development of neuropathic pain (Knerlich-Lukoschus et al., 2011), which is a condition that affects many SCI patients (Jutzeler et al., 2016). However, the clinical examination of pinprick sensation per se does not allow assessment of the directionality of this relationship. Therefore, it remains unclear whether the structural changes within the PAG are a result of primary loss of first order neurons or due to secondary processes affecting second order neurons, or both. Understanding the complex interaction of nociceptive (mal-) processing requires multimodal studies including direct readouts of nociceptive information flow (i.e. contact heat evoked potentials (Haefeli et al., 2013)) and structural as well as functional MRI. These multimodal biomarkers could then inform biophysical models of changes of states, which could be used to assess the directionality and extent of how the coupling between spinal cord, brainstem and brain is influenced by pain .

Iron content in brainstem pathways and nuclei
Interestingly, we did not observe iron accumulation caused by myelin breakdown (Hametner et al., 2013) causing detrimental inflammation in the CNS (Felix et al., 2012;Kroner et al., 2014). This may be due to rather small effects of iron accumulation in supraspinal regions compared to the spinal cord, where more iron is released by breakdown of haemoglobin after haemorrhage. Longitudinal assessment of iron accumulation after acute SCI will shed more light into these mechanisms with greater sensitivity to subtle effects.

Relationship between severity of injury and neurodegeneration
We did not observe either lesion level (tetra-/paraplegia) or severity (completeness) dependent neurodegeneration in brainstem sub-structures. This is in accordance with similar trajectories of atrophy above the level of injury between tetraplegic and paraplegic patients within the first year after injury within the spinal cord and brain (Freund et al., 2013). However, in chronic SCI lesion height and severity determined the magnitude of cord atrophy (Jutzeler et al., 2016). It therefore remains speculative why a cervical injury with a greater impact on the structural integrity of a higher number of fibres and neurons than a comparable thoracic lesion would not lead to more neurodegeneration within the brainstem. One potential explanation could be a flooring effect of atrophic progression due to the long-standing injuries. Moreover, given that the effect of trauma per se is strong, inducing linear and non-linear changes across the entire neuroaxis (Freund et al., 2013;Grabher et al., 2015), such changes, as well as differences in other factors (e.g. treatments, time spent in rehabilitation) and their complex interactions, may have concealed effects between patient sub-groups (i.e. tetra-/paraplegia, completeness). Finally, the cross-sectional design restricts conclusions to a single time point and fully characterizing these dynamic processes requires longitudinal studies. However, the results of the present study motivate the design of longitudinal studies with sophisticated neuroimaging protocols to establish these clinicopathological relationships.
Of note, quantitative measures of magnetisation transfer saturation and longitudinal relaxation rate provide information about macromolecular content within the neuronal tissue and are thus indirect measures sensitive to myelin. Post-mortem validation studies have shown the high association between MT-based measures and myelin as its main contributor (Schmierer et al., 2004;Turati et al., 2015). Although macromolecular content is the main component in R1, other contributors are water and iron content Rooney et al., 2007). Finally, both quantitative readouts (MT and R1) are sensitive but not specific to a single underlying pathological mechanism. Thus they complement each other by providing insights into different disease processes (e.g. atrophy, change in myelin and water content, iron deposition) in brainstem sub-structures (Felix et al., 2012;Freund et al., 2013;Grabher et al., 2015;Kroner et al., 2014).

Limitations
We note the following limitations of this study. Patients were on average 7.8 years older than controls. Furthermore, age was significantly different between groups for pre-upgrade, but not for postupgrade data. We therefore included age as a nuisance variable in all Fig. 4. Correlation between microstructural integrity within the brainstem nuclei and neurological and functional impairment. Regression models from extracted peak-voxel within the significant cluster are shown for illustrative purposes only (not adjusted for age, scanner, and total intracranial volume). P. Grabher et al. NeuroImage: Clinical 15 (2017) 494-501 statistical models to account for linear age-related effects. The scanner was upgraded during the study period (from Verio to Skyra fit ). Data from both patients and controls were acquired on both systems (the same ratio in each group) to minimize potential confounding effects due to the scanner upgrade. In addition, we used the scanner upgrade as nuisance variable in all GLMs. Due to the exploratory nature of this study, no adjustments for the number of contrasts tested were performed (Bender and Lange, 2001). The anatomical locations of brainstem ROIs (i.e. TPMs) and findings in SCI compared to healthy controls were carefully confirmed using a high-field MRI atlas (Naidich et al., 2009). The lack of specificity for pathological mechanisms influencing relaxation times and the small anatomical structures relative to the acquired resolution may conceal small effects and may be solved by high-resolution multiparametric MRI techniques . The latter will help to improve segmentation of brainstem substructures and increase specificity for brainstem pathways and nuclei. The cross-sectional study design allows us to only assess differences in MRI readouts between SCI and healthy controls, but not the underlying trajectories of structural changes. This design is also less sensitive, as higher between-subject variability may conceal weak effects (e.g. no observed differences between tetraplegic and paraplegic patients, no dependence on lesion level and severity). Of note, lateralized findings are often introduced by thresholding statistical parametric maps to account for multiple comparisons. To overcome this limitation, we aim to develop a longitudinal analysis pipeline to assess trajectories of structural change in brainstem pathways and nuclei in acute SCI.

Conclusion
Refined qMRI methods enable tracking of spatially specific neurodegeneration and structural reorganization in the brainstem following traumatic SCI. Next to measures of atrophy, which are rather insensitive and unspecific to the underlying pathology, we show myelin reductions across the brainstem. Therefore, these clinically relevant structural brainstem alterations, obtained with a qMRI protocol, could serve as neuroimaging biomarkers to monitor treatment efficacy and complement clinical assessments in clinical trials following SCI.

Funding
This study was funded by the Clinical Research Priority Program "NeuroRehab" of the University of Zurich and Wings for Life [WFL-CH-007/14]. The Wellcome Trust Centre for Neuroimaging is supported by core funding from the Wellcome Trust [091593/Z/10/Z]. Open access of this publication was supported by the Wellcome Trust.

Conflicts of interest
We declare no conflicts of interest.

Contributions
All authors were involved in the overall study design. PG analysed the data. CB and JA were involved in the methodological aspects of the study. PG and PF wrote the paper. All co-authors reviewed the paper.