Microstructural plasticity in nociceptive pathways after spinal cord injury

Objective To track the interplay between (micro-) structural changes along the trajectories of nociceptive pathways and its relation to the presence and intensity of neuropathic pain (NP) after spinal cord injury (SCI). Methods A quantitative neuroimaging approach employing a multiparametric mapping protocol was used, providing indirect measures of myelination (via contrasts such as magnetisation transfer (MT) saturation, longitudinal relaxation (R1)) and iron content (via effective transverse relaxation rate (R2*)) was used to track microstructural changes within nociceptive pathways. In order to characterise concurrent changes along the entire neuroaxis, a combined brain and spinal cord template embedded in the statistical parametric mapping framework was used. Multivariate source-based morphometry was performed to identify naturally grouped patterns of structural variation between individuals with and without NP after SCI. Results In individuals with NP, lower R1 and MT values are evident in the primary motor cortex and dorsolateral prefrontal cortex, while increases in R2* are evident in the cervical cord, periaqueductal grey (PAG), thalamus and anterior cingulate cortex when compared with pain-free individuals. Lower R1 values in the PAG and greater R2* values in the cervical cord are associated with NP intensity. Conclusions The degree of microstructural changes across ascending and descending nociceptive pathways is critically implicated in the maintenance of NP. Tracking maladaptive plasticity unravels the intimate relationships between neurodegenerative and compensatory processes in NP states and may facilitate patient monitoring during therapeutic trials related to pain and neuroregeneration.


INTRODUCTION
The underlying pathophysiology of neuropathic pain (NP) after spinal cord injury (1) is complex and involves alterations within 'bottom-up' nociceptive processing (ie, afferent integrity) and 'topdown' endogenous pain modulation. 1 Functional changes within descending modulatory pathways have been related to the presence of NP after spinal cord injury (SCI). 2 The descending pain modulatory network encompasses a cortical-subcortical-brainstem network involved in the modulation of afferent nociceptive information. 3 Key constituents involved in nociceptive information processing are the posterior insula, thalamus, periaqueductal grey (PAG), the anterior cingulate cortex (ACC) and dorsolateral prefrontal cortex (DLPFC) (for review, see Wiech 4 ). Within ascending nociceptive pathways, inhibition and facilitation of afferent input already occurs at the level of the dorsal horn where primary afferent fibres synapse onto projection neurons. 5 Modulatory changes within these regions can precipitate a pro-nociceptive state and potentially contribute to the emergence and chronification of NP. 6 Traumatic SCI triggers a cascade of traumainduced secondary neurodegenerative processes, involving demyelination and iron accumulation in the spinal cord and brain, 7 which can be tracked using quantitative MRI (qMRI). 8 To date, increases and decreases in brain and cord macrostructural (ie, volumetric) changes have been associated with the occurrence of NP. 9 Moreover, macrostructural tissue sparing (ie, ventral tissue bridges) could be related to the emergence and maintenance of NP after SCI. 10 However, microstructural correlates of these NP associated processes within areas undergoing atrophy and beyond are understudied. Based on recent studies illustrating that activity-dependent plasticity can translate into changes in myelin architecture, we hypothesised that such pathophysiological processes will be reflected in microstructural changes to myelin and iron content along nociceptive pathways.
We used a multi-parametric mapping (MPM) protocol 8 that provides contrasts that are sensitive to myelination (via magnetisation transfer (MT) saturation, longitudinal relaxation (R1) 11 12 ), and sensitive to iron content in tissue as well as in blood (using the effective transverse relaxation rate (R2*) 13 ) to track the complex relationship between structural and metabolic changes along the trajectories of nociceptive pathways and its relation to the presence and intensity of NP. To characterise simultaneously NP related changes across the neuraxis, we used a combined brain and spinal cord template 14 embedded in the statistical parametric mapping (SPM) framework. Finally, we applied multivariate source-based morphometry (SBM), which estimates interrelationships among voxels across the neuraxis to identify naturally grouped patterns of structural variation between groups. 15 The multivariate (SBM) tests for group effects are, in principle, much more sensitive than the equivalent mass univariate (voxel-based morphometry (VBM)) tests one would obtain from analysing the volumetric and microstructural images directly.

Spinal cord injury
This is because multivariate analyses do not try to assign a significance to each voxel or region-they test for distributed effects that covary among individuals. 16 Due to SBMs sensitivity to minute changes, this approach was chosen to detect microstructural changes. To dissociate NP-associated changes from traumainduced changes we first compared MRI indices between SCI patients with and without NP and then between SCI patients and healthy controls.

MATERIALS AND METHODS Participants
Thirty chronic traumatic SCI patients (13 with and 17 without NP (table 1)) and 23 healthy controls participated in the study (online supplemental S1). These participants were also enrolled in a longitudinal study. [17][18][19] Eligible individuals with a traumatic SCI and healthy controls were older than 18 years, and had no history of head and brain lesions, no pre-existing neurological, mental or medical disorders affecting outcome, and no contraindications to MRI.

CLINICAL ASSESSMENT
The neurological examination was performed according to the International Standards for Neurological Classification of Spinal Cord Injury. 20 The pain intensity was quantified using an 11-point numeric rating scale with '0' indicating no pain and '10' indicating the worst imaginable pain. NP was defined according to the standard taxonomy related to SCI, 21 which required its presence in an area of sensory deficit (at or below the lesion level) and other causes like musculoskeletal pain had to be ruled out during the clinical assessment. We used ANOVA tests to assess group differences for the demographics (ie, lesion level, American Spinal Cord Injury Association (ASIA) Injury Severity at baseline (AIS) grade).

DATA ACQUISITION
Structural whole-brain data, including the cervical cord down to vertebra C5, were acquired on a 3T Magnetom-Verio MRI scanner (Siemens Healthcare, Erlangen, Germany) equipped with a 16-channel radiofrequency (RF) receive head and neck coil and RF body transmit coil. The scanner was upgraded during the study period (from Verio to Skyra fit ) which resulted in 15 controls and 20 patients being scanned on the Verio and 8 controls and 10 patients being scanned on the Skyra system.
A whole-brain multi-parameter mapping (MPM) qMRI protocol, 22 using a multi-echo 3D FLASH (fast low-angle shot) sequence, was performed. The following parameters were used: FoV of 240×256×176 mm 3

DATA ANALYSIS
We used SPM12 spatial routines on the T1w images and MPM maps. As a novelty, we used a brain plus spinal cord (BSC) template 23 covering the brain and the upper cervical spinal cord in the unified segmentation approach, 24 which allowed us to assess brain and cervical cord changes within the same statistical framework. This BSC template was specified as a tissue probability map (ie, spatial prior) in the unified segmentation step. This procedure assumes every voxel to be drawn from an unknown mixture of 12 Gaussians, which were grouped into seven distinct tissue classes: grey matter (GM, using one Gaussian), white matter (WM, using one Gaussian) and cerebrospinal fluid (using one Gaussian), bone and air (using three Gaussians), soft tissue (using three Gaussians), non-neural tissues (using three Gaussians) and air/background (using three Gaussians). 24 For volumetric analysis, we applied SPM12's unified segmentation, including the BSC template, to each subject's MPRAGE image (ie, T1w). T1w images were non-linearly transformed into standard MNI (Montreal Neurological Institute) space using diffeomorphic group-wise registration (Dartel) implemented in SPM12. 25 This defined the space for subsequent modelling steps. Finally, GM and WM probability maps were smoothed with an isotropic Gaussian kernel of 5 mm full width at half maximum.
For microstructural analysis, the acquired T1w, PDw and MTw FLASH echoes from the MPM protocol, were first averaged to increase signal to noise ratio. Unified segmentation based correction of R1 brain maps for RF transmit field inhomogeneities (UNICORT) was then used for correcting RF transmit field inhomogeneity and to calculate quantitative maps of MT saturation, R1 and R2*. 8 MT maps were segmented using the unified segmentation approach that included the BSC template. 14 Segmented MT maps were then transformed nonlinearly to standard MNI space using Dartel, but without scaling by the Jacobian determinants (ie, no 'modulation' of parameter values). Additionally, all MPM maps (MT, R1 and R2*) were warped to the MNI space by using the participant specific flow fields from the MT maps and finally smoothed using a previously established tissue-weighted-smoothing procedure with a kernel of 5 mm, in order to preserve quantitative values within the GM and WM tissue classes of each MPM map. Subsequent modelling and analysis were performed for smoothed, normalised morphometric and microstructural parameter images (for WM and GM) across the brain and cervical cord, simultaneously. Each of the morphometric and microstructural parameter images were then individually analysed using SBM. 16 SBM 16 is a multivariate analysis approach, implemented in the SBM toolbox within Group ICA Toolbox (GIFT) (http:// mialab. mrn. org) that estimates covarying networks. SBM can be considered as a multivariate extension of VBM which accounts for spatial dependencies among different regions and increases sensitivity to effects distributed across the neuraxis. SBM applies spatial independent component analysis 26 (organised as subjectsby-voxels) to produce maximally independent components (components-by-voxels) and their associated mixing-matrices (subjects-by-components). The number of components (k) for subsequent ICA were estimated using a principled information theoretic approach (rather than arbitrarily selecting the number of components 26 ). ICA was then applied to the T1w WM images of 23 HC, 12 NP and 17 without NP, which were arranged into one 53-row (subject-by-T1w-WMvoxels) data matrix. This data matrix was then decomposed into a mixing matrix (subjects by components) and a component matrix (components by voxels). The mixing matrix expresses the relationship between 53 subjects and k components. The rows of the mixing matrix indicate the contribution of each k component to a given subject, whereas the columns indicate how each component contributes to the 53 subjects. The component matrix on the other hand expresses the relationship between the k components and brain/cord voxels. The rows of the component matrix indicate how one component contributes to different voxels, whereas the columns of the component matrix indicate how one voxel contributes to each of the components.
The mixing matrix was used for subsequent statistical analysis, where every column of the mixing matrix quantifies the contribution of each component to the 53 subjects. An ANOVA was applied to each column to assess which components showed a significant group difference, adjusted for covariates of no interest (ie, age, gender, total intracranial volume (TIV) and scanner). A false discovery rate 27 controlling procedure with q*=0.05 was used to assess which components were statistically significant. The significant component that survived FDR correction was thresholded at a value of |Z|>3 as described in Xu et al 28 and was scaled to unit SD (SBM Z map). A significant component comprises clusters of voxels-with positive and negative values-reflecting group effects in different regions of brain and spinal cord that are interrelated to each other.
SBM was repeated for each of the volumetric and microstructural parameter maps. The number of components estimated were 8, 8, 7, 7, 6, 3, 5, 5 for GM and WM maps from the MRPAGE, R1, R2* and MT maps, respectively. Components with loading scores that differed significantly between either controls and patients or NP patients versus pain-free patients were identified. No significant components were found for MT maps. The significant components were rendered on the MNInormalised BSC template for appropriate contrasts (controls>patients, patients>controls, NP patients>pain free, pain-free>NP patients) (figures 1 and 2). For illustrative purposes, the GM and WM are represented in the same colour for each modality (eg, volume, R1, R2*) of the significant components.
Finally, in order to explore the association between brain and spinal cord changes and pain intensity in patients with NP, we used SPM's multiple linear regression models adjusting for potentially confounding effects of age, gender, time since injury, TIV and scanner upgrade. The explanatory variables in these analyses were pain intensities (range from 0 to 10), while the response variables were the volumetric and microstructural parameter maps above.

Patients' characteristics and clinical outcomes
The demographics of the patients is shown in table 1. Out of 30 chronic traumatic SCI patients (3.6±6.8 years after the injury), 13 reported NP (mean age: 42.9 years, SD=18.2, range=21-73) and 17 were pain-free (mean age: 46.1 years; SD=14.8, range=19-72). All patients underwent state of the art neurorehabilitation in a specialised SCI centre immediately after their injury. The healthy controls had a mean age of 36.9 years (SD=11.8, range=24.0-66.0). All three groups were comparable with respect to age (ANOVA, p=0.1516) and patient groups did not differ with regard to lesion level (p=0.28) or AIS grade (p=0.18). The mean below level NP intensity was 6.1 (SD=1.9).

Pain-associated changes in SCI patients with NP compared with pain-free SCI patients
In NP patients, atrophic changes associated with NP were observed in the cerebellum, middle temporal gyrus and occipital gyrus, while volume increases were observed in the posterior insula, superior temporal, occipital and middle frontal gyrus when compared with pain-free patients (figure 1A and online supplemental S2). Analysis of the R1 component showed signal decreases in primary motor cortex and DLPFC in NP patients when compared with pain-free patients (figure 2A and online supplemental S3). In NP patients, the R2* component showed increased signal in the cervical cord, thalamus, ACC, PAG and para-hippocampal gyrus, while signal was decreased in the lentiform nucleus and substantia nigra when compared with painfree patients ( figure 2A and online supplemental S4). Note, that trauma-related differences were predominantly visible in the medial thalamus whereas pain-related changes were located in the lateral thalamus.

Trauma-associated changes in SCI patients compared with healthy controls
We first confirmed previous reports showing trauma-induced macrostructural and microstructural changes across the motor, Spinal cord injury sensory and limbic system. Specifically, macroscopic changes (decreased volume) were observed in the cervical cord, cerebellum, thalamus, posterior cingulate, lingual gyrus, precuneus, superior and medial temporal gyrus and para-hippocampal gyrus (figures 1B and 2) when all patients (NP+pain free) were compared with controls. In these atrophied areas and beyond,

DISCUSSION
This study shows both decreased and increased myelin-sensitive and iron-sensitive content differences associated specifically with NP along the trajectory of ascending and descending nociceptive pathways after SCI. Most group effects are located within key structures involved in endogenous pain modulation. 4 As SCI provides an ideal model to differentially study trauma-related and pain-related effects-following a defined lesion to the nervous system-these findings speak to specific pathophysiological mechanisms related to the development and maintenance of chronic NP after SCI.
We first dissociated trauma-induced changes from NP associated (micro-)structural changes. In agreement with previous reports, trauma-induced atrophic changes, as well as myelin and iron changes, were evident across the sensorimotor and limbic system (for review, see Freund et al 7 ). Accompanying trauma-induced changes, volumetric increases and decreases in the cervical cord, thalamus, primary somatosensory cortex, ACC and DLPFC have been also associated with NP following SCI. 9 The neuronal substrates underlying injury-induced volumetric alterations remain incompletely understood. Pathophysiologically, chronic ectopic electrical activity in residual spinothalamic tract axons, 29 upregulation of sodium channels in nociceptive neurons, 30 and increases in pro-inflammatory cytokines 31 are potentialmechanism involved. As volumetric changes Figure 2 Source-based morphometry renderings showing trauma induced macrostructural and microstructural changes across the motor, sensory and limbic system (|Z|=3) in the brain and spinal cord. (A) Trauma induced volumetric changes (ie, all patients (neuropathic pain (NP)+pain-free) were compared with controls) the cervical cord, thalamus, posterior cingulate, lingual gyrus, precuneus, superior and middle temporal gyrus and para-hippocampal gyrus showed significant atrophy (shown in blue). (B) For patients with and without NP when compared with healthy controls, myelin-sensitive longitudinal relaxation (R1) component (shown in green) was decreased in the sensorimotor cortices, middle frontal gyrus, precuneus and inferior parietal lobule, whereas iron-sensitive effective transverse relaxation rate (R2*) component (shown in pink) showed increases in the cervical cord, lentiform nucleus, cerebellum, posterior cingulate, middle temporal gyrus, right cerebellum and fusiform gyrus. MPRAGE, magnetisation-prepared rapid acquisition gradient echo; T1w, T1 weighted.

Spinal cord injury
can represent a myriad of physiologically and pathophysiological changes, we sought to reveal the underlying microstructural substrates of such pain-associated volumetric changes by revealing changes in myelin and iron content along the trajectories of ascending and descending nociceptive pathways. Specifically, R2* increases were detectable in the spinal cord, PAG, thalamus, ACC and para-hippocampal gyrus, suggesting iron accumulation due to neurodegenerative processes 13 in regions that have previously been implicated in NP after SCI. 9 18 32 It remains speculative which pathophysiological processes drive changes in R2* contrasts as a range of physiological and pathophysiological processes can contribute to the gradient echo signal decay (and thus R2*). Remyelination, during recovery and reorganisation processes may impact myelin architecture and thereby R2*. 13 Inflammatory processes leading to microglial activation and neurodegenerative processes may result in changes in R2*, but also to changes in metabolic demand affecting the concentration of deoxyhaemoglobin in certain brain areas also influences R2*. 33 However, although there is a complex interplay between myelin, iron accumulation and R2* values in brain pathology, remote spinal cord and brain changes following SCI are likely related to anterograde and retrograde degeneration and subsequent changes in diaschitic regions.

Figure 3
Statistical parametric maps testing (ie, linear regression) for an association between longitudinal relaxation (R1) and effective transverse relaxation rate (R2*) signal changes and pain intensity in the periaqueduct and medulla oblongata (t values uncorrected p<0.005, shown for descriptive purposes, masked by the union of the regions of interest). PAG, periaqueductal grey.
Communication between the prefrontal cortical areas and the PAG has been recently shown in functional MRI studies. 34 Here, we provide evidence for further downstream effects. Notably, microstructural alterations were detected within the spinal cord and also the PAG. These effects were related to pain intensity; furnishing an important predictive validity to this imaging phenotype. The PAG is known to serve as a link between the forebrain and the lower brainstem and is pivotal for the descending modulation of pain. It receives input from the frontal lobe, amygdala, hypothalamus and ACC and engages the rostral ventral medulla which in turn controls nociceptive processing in the spinal dorsal horn. 35 Direct stimulation of the PAG has successfully been trialled as a treatment option for NP including pain following SCI. 36 Thus, it is anticipated that molecular changes in the PAG could compromise descending inhibition of nociceptive processing and subsequently increase the spinal gain of incoming information. Our study is the first to provide evidence for NP-related changes in PAG microstructure in humans, underscoring the clinical relevance of these changes for NP following SCI. In conjunction with microstructural alteration in the injured spinal cord, these findings are compatible with the notion of a dysfunctional spinal-bulbo-spinal loop, a key circuit of the descending pain modulatory system. 3 Additional evidence for an impaired descending pain modulatory system can be inferred from microstructural changes in other key constituents of this network, namely the ACC (ie, increased iron accumulation) and DLPFC circuitry 4 (ie, reduced myelin content). Previous studies had found structural and functional NP associations in the DLPFC 37 in those with NP. For example, reduced GM volume and hypometabolisms in the left DLPFC were observed in SCI patients suffering from NP when compared with healthy controls. 38 Here, we not only show that this reduction in volume is linked to a decrease in myelin-sensitive R1, but also that it is pain specific-as this reduction was not observed in pain-free patients. As the DLPFC is both involved in top-down inhibition and facilitation of pain, the decreases in myelin content reported herein align well with the concept of pathological nociceptive gain control in NP after SCI. 2 The thalamus is considered to play an important role in the pathophysiology of central NP after SCI. 32 After a spinal lesion ascending somatosensory pathways are damaged, which may result in deafferentation of rostral relay structures such as the thalamus. Intriguingly, our finding of R2* signal changes which may reflect iron accumulation in the thalamus in SCI patients with NP was confined to lateral thalamic nuclei, while volumetric changes in the medial thalamus were trauma-related (ie, found in SCI with and without NP). These findings offer an essential line of evidence that structural plasticity within the thalamus is specifically linked to the pathophysiology of NP after SCI.
Finally, we observed an unexpected decrease of iron-sensitive R2* in the basal ganglia. The role of the basal ganglia in chronic pain conditions is still poorly understood, but its vast anatomical connections to a multitude of brain areas (including the thalamus) and cervical cord make it plausible that the basal ganglia could play an important role in aberrant nociceptive bottom-up and top-down signalling (for review, see Borsook et al 39 ). Our findings warrant further studies into the role of basal ganglia pathways in central NP. At present, the implications of these findings remain highly speculative.
From a technical perspective, the multivariate (SBM) tests for group effects described above are, in principle, much more sensitive than the equivalent mass univariate (VBM) tests one would obtain from analysing the volumetric and microstructural images directly. This is because multivariate analyses do not try to assign a significance to each voxel or region-they test for distributed effects that covary among individuals. 16 Crucially, this multivariate analysis would not have been possible without combining the brain and spinal cord within the same analysis. This speaks to the potential importance of the combined brain and spine template used to spatially normalise our data. This template is available through open access (http://www. fil. ion. ucl. ac. uk/ spm/ toolbox/ TPM/) for related studies that test for distributed effects throughout the neural axis.

LIMITATIONS
Our study had some limitations. The cross-sectional nature of the study restricts conclusions to a single time point and thus the temporal evolution of the above-described microstructural changes remains uncertain. Despite the histological evidence that MT, R1 and R2* markers correspond to their biochemical counterparts 13 , they are indirect contrasts of myelin and iron content and any interpretation should take this into consideration.
It should be noted that concurrently observed changes in R2*, R1 and MT may not be apparent in all microstructural changes. First, the sensitivity or signal-to-noise ratio in these different quantitative measures varies significantly. 40 Thus, differences in sensitivity may render these measures complementary, that is, certain changes may only be visible in one of the metrics. Second, the different metrics have a distinct specificity and sensitivity to underlying microstructural changes. 22 For example, R2* is exquisitely sensitive to changes in local susceptibility and concentration of paramagnetic compounds such as iron. Callaghan and colleagues 41 have demonstrated that R1 can be estimated from R2* and MT measurements by a general linear relaxometry model, which demonstrates that R2* and MT changes may even cancel each other and result in zero change in R1. Thus, partial contributions of unexplored physiological/cellular processes occurring after SCI cannot be excluded. Moreover, current standardised neurological tests cannot account for unobserved latent lifestyle or genetic factors which might be different between SCI patients and controls, a-priori. To mitigate any potential effect of the scanner upgrade on our results, we ensured that the same number of patients and controls were measured before and after upgrade, allowing us to account for the upgrade effect (common to both cohorts) and modelled out any linear effect of these covariates of interest. Moreover, Leutritz et al 40 demonstrated that the systematic bias in R1, MT and R2* measurements between a Skyra and Verio Siemens scanner setup is <4%. Thus, we are confident that the effects seen in this study are related to pathophysiological changes rather than technical sources. Finally, sex was not balanced across groups, with most of the participants being men. However, this is representative of the general population of SCI patients, in which the male to female patients' ratio is roughly 4:1 and we also adjusted our models for this covariate of no interest.

CONCLUSION
This study evinces the microstructural signature of NP, affecting key constituents of the ascending and descending nociceptive pathways-its magnitude being directly linked to NP intensity. The complex interplay between myelin and iron changes in areas related to sensory and affective processing highlights maladaptive plastic processes likely involved in the maintenance of NP. Beyond unravelling the intimate pathophysiology of NP, tracking microstructural plasticity may facilitate patient monitoring during clinical trials for NP.