Multiple sclerosis lesions affect intrinsic functional connectivity of the spinal cord

Multiple sclerosis lesions affect intrinsic functional connectivity of the spinal cord Abstract Patients with multiple sclerosis present with focal lesions throughout the spinal cord. There is a clinical need for non-invasive measurements of spinal cord activity and functional organization in multiple sclerosis, given the cord’s critical role in the disease. Recent reports of spontaneous blood oxygenation level-dependent fluctuations in the spinal cord using functional MRI suggest that, like the brain, cord activity at rest is organized into distinct, synchronized functional networks among grey matter regions, likely related to motor and sensory systems. Previous studies looking at stimulus-evoked activity in the spinal cord of patients with multiple sclerosis have demonstrated increased levels of activation as well as a more bilateral distribution of activity compared to controls. Functional connectivity studies of brain networks in multiple sclerosis have revealed widespread alterations, which may take on a dynamic trajectory over the course of the disease, with compensatory increases in connectivity followed by decreases associated with structural damage. We build upon this literature by examining functional connectivity in the spinal cord of patients with multiple sclerosis. Using ultra-high field 7 T imaging along with processing strategies for robust spinal cord functional MRI and lesion identification, the present study assessed functional connectivity within cervical cord grey matter of patients with relapsing-remitting multiple sclerosis (n = 22) compared to a large sample of healthy controls (n = 56). Patient anatomical images were rated for lesions by three independent raters, with consensus ratings revealing 19 of 22 patients presented with lesions somewhere in the imaged volume. Linear mixed models were used to assess effects of lesion location on functional connectivity. Analysis in control subjects demonstrated a robust pattern of connectivity among ventral grey matter regions as well as a distinct network among dorsal regions. A gender effect was also observed in controls whereby females demonstrated higher ventral network connectivity. Wilcoxon rank-sum tests detected no differences in average connectivity or power of low frequency fluctuations in patients compared to controls. The presence of lesions was, however, associated with local alterations in connectivity with differential effects depending on columnar location. The patient results suggest that spinal cord functional networks are generally intact in relapsing-remitting multiple sclerosis but that lesions are associated with focal abnormalities in intrinsic connectivity. These findings are discussed in light of the current literature on spinal cord functional MRI and the potential neurological underpinnings. multiple sclerosis, spinal cord, lesions, functional connectivity, resting state functional magnetic resonance imaging Introduction Clinical background The spinal cord is a critical structure in health and disease. Damage to the cord can lead to dysfunction of motor and/or sensory systems, as occurs in cases of trauma involving spinal cord injury, as well as in various musculoskeletal and neurodegenerative diseases such as degenerative disc disease and multiple sclerosis. Spinal cord dysfunction may also be implicated in the absence of gross structural damage, such as in chronic pain syndrome. Given the central role of spinal cord dysfunction in these patient populations, there is an important (but largely unmet) clinical need for non-invasive measurements of cord activity and functional organization (Wheeler-Kingshott et al., 2014). Functional connectivity in the spinal cord We recently reported the presence of resting state blood oxygenation level-dependent (BOLD) correlations in the human spinal cord using ultra-high field MRI, as well as evidence that these effects are reproducible within and across human subjects (Barry et al., 2014, 2016). Such observations have since been corroborated by several research groups using alternative acquisition, processing, and statistical methods (Kong et al., 2014; Liu et al., 2016; San Emeterio Nateras et al., 2016; Eippert et al., 2017b; Harita and Stroman, 2017; Weber et al., 2018), as well as in non-human primate and rat spinal cord (Chen et al., 2015; Wu et al., 2018). These studies have demonstrated that spontaneous BOLD signal fluctuations in the spinal cord are primarily characterized by ventral motor and dorsal sensory networks. The organization and spatial distribution of this activity at rest is also consistent with task-based functional MRI studies (Maieron et al., 2007; Cadotte et al., 2012; Stroman et al., 2012). Eippert et al. (2017b) looked at resting state correlations in cervical cord grey matter at 3 T using multiple variations on their preprocessing and analysis pipelines, finding that the ventral-ventral and dorsal-dorsal horn networks were robust to variations in processing procedures, successfully replicating our findings at 7 T. However, it is important to examine how spontaneous spinal cord activity and functional connectivity is altered in disease states, particularly diseases that result in cord dysfunction. Overall goals In the current study, we (i) examined resting state BOLD fluctuations in the cervical spinal cord grey matter in patients with relapsing-remitting multiple sclerosis compared to a large sample of healthy controls; (ii) assessed the relationship of spinal cord connectivity measures with standard clinical measures of disability; and (iii) assessed the effects of lesions on region-to-region connectivity both locally and distant in relation to ascending/descending white matter tracts. The present study provides the first assessment of functional MRI-derived functional connectivity in the cord of a diseased population. Structural and functional organization of the spinal cord The spinal cord is responsible for the transmission and conduction of electrical potentials to and from the brain, the primary interface of the central and peripheral nervous systems. It is well organized and is composed of grey and white matter regions. The spinal cord grey matter is centrally located surrounded by white matter and is characterized by its ‘butterfly’ shape. Spinal cord white matter contains dense bundles of myelinated axons running in the rostrocaudal direction, separated into the dorsal, lateral, and ventral columns. Spinal cord grey matter can be divided into functional subunits corresponding to motor and sensory inputs, including anterior (ventral) and posterior (dorsal) horns (Cohen-Adad and Wheeler-Kingshott, 2014). Motor network Descending spinal pathways originating primarily in motor cortex are responsible for voluntary motor control. Projections from contralateral cortical sites decussate (cross hemispheres) in the medulla of the brainstem and travel down the cord along ventral and lateral corticospinal tracts. Corticospinal tract axons synapse onto neurons in the ventral grey matter of the spinal cord, from which efferent fibres exit the ventral roots of the cord and innervate specific muscle groups (Purves et al., 2012; Kandel et al., 2013). Sensory network Ascending spinal pathways carry somatosensory information up to the brain from a myriad of peripheral afferents, including tactile, temperature, and pain receptors in the skin as well as from muscles, joints, and visceral organs. Mechanoreceptors and proprioceptors project onto the dorsal column medial lemniscus (DCML) pathway and are responsible for encoding information about fine touch, vibration, pressure, and position. These receptors send their first-order projections from ipsilateral dorsal root ganglia (location of sensory receptor cell bodies) to the spinal cord where the axons travel along ipsilateral dorsal column tracts and ultimately synapse onto grey matter neurons in the medulla. Alternatively, nociceptors and thermoreceptors, which encode information about pain, temperature, and coarse touch, project onto the spinothalamic pathway, which involves anterolateral white matter tracts. For spinothalamic afferents, first-order axons from ipsilateral dorsal root ganglia enter the spinal cord and then synapse onto grey matter neurons in the ipsilateral dorsal horn of the cord. Interneurons and central pattern generators Local connections exist within spinal cord segments that may demonstrate distinct, synchronized activity. The existence of central pattern generators (CPGs) in the spinal cord was observed over 100 years ago (Brown, 1911; McCrea and Rybak, 2008). In the spinal cord, CPGs refer to neural computations of locomotor actions and reflexes, which are performed within the spinal cord and without specific input from higher-order brain areas. These computations are mediated by interneurons that connect local grey matter regions, allowing, for instance, sensory input to locally influence ventral motor neuron output. Challenges in spinal cord functional MRI While functional MRI studies of the brain have increased dramatically over the past 25 years, the detection of robust BOLD signal in the spinal cord has proved challenging (Giove et al., 2004; Stroman et al., 2014). The cervical spinal cord is ∼1–1.5 cm in diameter and grey matter butterfly horns only a couple of millimetres across, and the typical spatial resolution of 2–3 mm used in brain functional MRI is inadequate for spinal cord functional MRI. Higher in-plane resolution is needed to accurately distinguish spatial properties of BOLD signals across the spinal cord and to reduce partial volume effects from nearby tissue such as white matter and CSF. Additionally, susceptibility distortions arising from single-shot echo-planar imaging (EPI), used in brain functional MRI sequences are amplified in spinal cord due to the close proximity to surrounding bones and intervertebral discs and the oesophagus and lungs. A further complication in measuring the spinal cord BOLD signal is an increased influence of physiological noise compared to the brain (Eippert et al., 2017a). These noise sources include: cord motion, respiration and swallowing, transient B0 inhomogeneity, flow and CSF pulsation. We have developed improved acquisition and processing strategies for robust spinal cord functional MRI (Barry et al., 2011, 2013, 2014, 2016, 2018; Conrad et al., 2017), which are used in this study and outlined in the ‘Materials and methods’ section. Functional connectivity of sensorimotor networks in multiple sclerosis Given the spinal cord’s high degree of axonal connections to sensorimotor cortices, as well as reports of synchronous task-based BOLD activity in spinal cord with these brain regions, it is likely that spinal cord resting state BOLD activity is associated with sensorimotor network activity. Studies looking at functional connectivity in multiple sclerosis in the brain have demonstrated a number of abnormalities compared to healthy controls, with reports of both increased and decreased correlation among networks. Further evidence suggests that these abnormalities may be dynamic over time, with compensatory increases in connectivity occurring early in the disease and decreased connectivity with the accumulation of structural damage (Roosendaal et al., 2010; Faivre et al., 2012). It is possible that analogous changes occur in the sensorimotor network with improvement of motor deficits in response to treatment, for example. Such improvements may manifest in and/or involve spinal cord circuits. While more work is needed, the growing body of evidence supporting widespread network alterations in multiple sclerosis suggests resting state functional MRI measures may provide useful indicators of disease pathology and progression (Rocca et al., 2012b, 2015; Filippi and Rocca, 2013; Sacco et al., 2013). Complementary to studies of the brain, the investigation of intrinsic BOLD fluctuations and functional connectivity in the spinal cord may thus provide a more comprehensive understanding of nervous system dysfunction in multiple sclerosis. Furthermore, measures of functional connectivity in the cord may provide surrogate markers of sensorimotor deficits and aid in disease monitoring or prognosis. Functional imaging of the spinal cord in multiple sclerosis While we present new findings in spinal cord resting state functional MRI in patients with multiple sclerosis, several previous studies have looked at task-based functional activation of the cervical cord in multiple sclerosis in response to tactile and proprioceptive stimulation of the hand or wrist (Agosta et al., 2008a, b, 2009; Valsasina et al., 2010, 2012). This body of work has shown that patients with multiple sclerosis demonstrate increased activation of the cervical cord compared to control subjects, as well as a more bilateral response in the context of ipsilateral stimulation. We may expect over-recruitment of bilateral networks in a task-paradigm to manifest as increased connectivity in a resting state, and therefore predicted that patients with multiple sclerosis would demonstrate an overall increase in connectivity between bilateral region pairs in comparison to controls. Importantly, abnormalities detected via spinal cord functional MRI may relate to clinical measures in multiple sclerosis. One study looking at the relationship of fatigue in patients with multiple sclerosis with cord activity in response to palm stimulation found increased activation in patients with low fatigue scores as compared to controls (Rocca et al., 2012a). Interestingly, a second group of patients with high fatigue scores showed reduced activation compared to those with low fatigue, and instead demonstrated a more diffuse pattern of activity compared to controls. Similar to findings in the brain in multiple sclerosis, it is possible that abnormal functional activity in the cord may follow a dynamic trajectory of early compensation followed by decreased recruitment in later stages concomitant with increasing disability and structural damage. In addition to assessing global abnormalities in spinal cord functional connectivity in multiple sclerosis, a primary aim of the present study was to assess the effects of spinal cord lesion location in multiple sclerosis on functional connectivity among grey matter regions of the cervical cord. Dula et al. (2016) recently showed that 50% more lesions were detectable with 7 T imaging compared to 3 T using similar scan times in the same group of patients. The present study takes advantage of this increased sensitivity to lesions through the use of 7 T imaging. In regards to prior work on spinal cord connectivity and sites of local damage such as with multiple sclerosis lesions, a study by Chen et al. (2015) looked at functional connectivity in the cervical cord of anaesthetized squirrel monkeys before and after introducing a unilateral dorsal lesion in the context of a spinal cord injury model. Connectivity was reduced within the slice of the lesion as well as slices below the lesion involving ipsilateral ventral and dorsal horns, whereas connectivity above the lesion and in contralateral horns was largely unaffected. These results demonstrate that focal lesions in the cord may have both local (intrasegmental) as well as distal (intersegmental) effects on spontaneous fluctuations in spinal grey matter activity. Recovery of disrupted connectivity towards pre-injury levels was observed in a time period of 7–24 weeks after lesioning. This suggests that reorganization of spinal circuitry occurs in response to damage and that monitoring of resting state connectivity may serve as a useful clinical indicator of functional integrity. It is important to note that lesions in the spinal cord are less well studied than those in the brain, largely due to limitations in the spatial resolution of MRI at clinical field strengths. Some work has shown that specific lesion locations, such as within the corticospinal tracts and motor network in general, are more predictive of disease progression (Bodini et al., 2011), providing further motivation for assessments of spinal cord lesions and their impact on functional integrity. Materials and methods Subjects and image acquisition All subjects provided signed, informed consent prior to inclusion in the study in accordance with the local institutional review board. Patients with relapsing-remitting multiple sclerosis (n = 22) were recruited in collaboration with the Vanderbilt Multiple Sclerosis Center as part of a multi-modal 7 T MRI study. Healthy volunteers (n = 56) were also recruited and the first resting state run was used for analysis if two were acquired. Subject demographics for our patient and control groups are provided in Table 1. An expanded set of clinical variables and symptoms for each patient is provided in Supplementary Table 1. All imaging was performed on a Philips Achieva 7 T scanner using a cervical spine coil (Nova Medical Inc.) with quadrature transmit and a 16-channel receive array. Foam padding was used in positioning subjects on the coil such that the cervical spine was straight but comfortable. Subjects were instructed to lie still throughout the scan. Physiological monitoring was performed for all subjects and involved recording from a pulse oximeter placed on a finger and respiration bellow placed on the abdomen. Table 1 Subject demographics Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  EDSS = Expanded Disability Status Scale; TSNR = temporal signal-to-noise ratio. Table 1 Subject demographics Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  EDSS = Expanded Disability Status Scale; TSNR = temporal signal-to-noise ratio. For each subject, a high resolution sagittal T2-weighted image was initially acquired and used to plan axial imaging such that the axial slice stack was centred on the C3-C4 vertebral disc and perpendicular to the spinal cord. Curvature of the cervical spine was sometimes not fully mitigated by padding and thus an orthogonal slice placement was not possible. In such cases, slice angulation was planned such that the middle slices were perpendicular to the cord at the C3-C4 vertebral level. Prior to functional imaging, a T2*-weighted anatomical volume was acquired with the same axial slice geometry as the functional scan and the following parameters: field of view = 160 × 160 mm2, 12 4-mm slices (centred on the C3/C4 junction), nominal voxel size = 0.6 × 0.6 × 4 mm3, interpolated voxel size = 0.31 × 0.31 × 4 mm3, repetition time = 303 ms, echo time = 8.2 ms, flip angle = 25°, sensitivity encoding (SENSE) (Pruessmann et al., 1999) reduction factor = 2.0 (right-left), number of acquisitions = 8, total acquisition time = 5 min and 22 s. The anatomical images provided high resolution detail of the cord, including grey matter, white matter, and CSF boundaries, and were used to delineate structures of interest for functional data processing. Additionally, the T2*-weighted anatomical images were used for lesion ratings in patients. For functional imaging, a T2*-weighted resting state functional MRI scan with 150 dynamics was acquired, using a 3D multishot gradient echo sequence with the following parameters: field of view = 160 × 160 mm2, 12 4-mm slices, voxel size = 0.91 × 0.91 × 4 mm3, repetition time = 17 ms, echo time = 8.0 ms, flip angle = 15°, echo train length (number of k-space lines acquired per shot) = 9, SENSE reduction factor = 1.56 (anterior-posterior), volume acquisition time = 3.34 s (278 ms/slice), number of volumes = 150 (after 10 ‘dummy’ scans), total scan time ∼9 min. Preprocessing was the same for each subject and followed a pipeline previously described in detail (Barry et al., 2016), involving motion correction, physiological signal de-noising, white matter and CSF signal regression, and bandpass filtering from 0.01–0.13 Hz (Cox, 1996; Glover et al., 2000). See the online Supplementary material for more information on these procedures. Region of interest correlation analysis The grey matter masks were divided into four quadrants (regions of interest) (Fig. 1) representing the ventral and dorsal horns, excluding the central grey matter commissure and central canal. The grey matter regions of interest were then morphologically eroded to remove the outermost layer of interpolated voxels, minimizing partial volume effects from white matter. The mean signal time course across all voxels within each region of interest was extracted. Pearson correlation coefficients were computed for each of the six region of interest pairs within each slice. The correlation coefficients were transformed to z-scores using a Fisher transformation: z = tanh−1(r)(dof − 3)1/2, where dof is the estimated degrees of freedom after correction for first-order autocorrelation (Rogers and Gore, 2008) and z represents the distance from zero correlation in units of standard deviation for each region of interest pair. Note that while we have explored the use of partial correlations in previous work (i.e. including as nuisance variables the time series from the regions of interest opposite to the pair being correlated), we chose to analyse these data with full correlations to simplify interpretation of the results. Figure 1 View largeDownload slide Imaging and resting state connectivity.Top left: Example slice placement as centred over C3-C4 vertebral disk, used as landmark for placement in all subjects. Top right: Example high resolution, T2*-weighted anatomical image and close-up; interp. = interpolated voxel size (displayed in figure). Middle right: Example functional MRI-weighted acquisition and schematic of functional imaging showing 150 volumes acquired over time (t), which was every 3.34 s. Functional images acquired using a 3D multishot gradient echo sequence, images shown after registration and interpolation to anatomical space. Middle left: Example seed-to-voxel correlation map (thresholded at Pearson’s r > 0.465, P < 10−7 uncorrected). An example seed voxel in left ventral grey matter (GM) shows correlation with contralateral ventral, ipsilateral dorsal, and contralateral dorsal grey matter. Bottom left: Spinal cord regions of interest defined during manual segmentation and preprocessing. CC = central commissure of grey matter (voxels excluded from analyses); LD = left dorsal grey matter horn; LV = left ventral grey matter horn; RD = right dorsal grey matter horn; RV = right ventral grey matter horn; NS = ‘not-spine’ mask (included all voxels outside of the spinal canal); WM = white matter. Bottom right: Representative functional MRI BOLD signal timecourses in a single slice from LV and RV, and LD and RD, respectively, after bandpass filter. Pearson’s r and associated Z-score for correlation between both region of interest pairs provided. Figure 1 View largeDownload slide Imaging and resting state connectivity.Top left: Example slice placement as centred over C3-C4 vertebral disk, used as landmark for placement in all subjects. Top right: Example high resolution, T2*-weighted anatomical image and close-up; interp. = interpolated voxel size (displayed in figure). Middle right: Example functional MRI-weighted acquisition and schematic of functional imaging showing 150 volumes acquired over time (t), which was every 3.34 s. Functional images acquired using a 3D multishot gradient echo sequence, images shown after registration and interpolation to anatomical space. Middle left: Example seed-to-voxel correlation map (thresholded at Pearson’s r > 0.465, P < 10−7 uncorrected). An example seed voxel in left ventral grey matter (GM) shows correlation with contralateral ventral, ipsilateral dorsal, and contralateral dorsal grey matter. Bottom left: Spinal cord regions of interest defined during manual segmentation and preprocessing. CC = central commissure of grey matter (voxels excluded from analyses); LD = left dorsal grey matter horn; LV = left ventral grey matter horn; RD = right dorsal grey matter horn; RV = right ventral grey matter horn; NS = ‘not-spine’ mask (included all voxels outside of the spinal canal); WM = white matter. Bottom right: Representative functional MRI BOLD signal timecourses in a single slice from LV and RV, and LD and RD, respectively, after bandpass filter. Pearson’s r and associated Z-score for correlation between both region of interest pairs provided. Lesion ratings Lesion ratings were independently conducted by three raters (B.N.C., S.M., and A.B.) who were blind to patient clinical information during the rating process. Control images of the spinal cord were also provided to the raters to determine the normal signal intensity of the axial T2*-weighted images, e.g. hyperintensities present in healthy white matter, in an effort to eliminate false-positive lesion identification in the patient data. For every slice (264 total slices across all patients), raters determined whether a lesion was present in each of the four white matter columns: ventral, dorsal, left lateral, and/or right lateral, yielding a total of 1056 possible lesion locations. Raters also had the option to not rate an entire slice if an artefact interfered with lesion demarcation. After independent ratings were completed, a consensus meeting was held in which slices were reassessed in the presence of all three raters. Discrepancies in the original ratings were reviewed and a final rating was agreed upon in each instance. Statistical analyses Group differences of connectivity values in patients compared to controls were assessed with two-tailed Wilcoxon rank-sum tests, performed using the mean connectivity value across all 12 slices for each subject. To assess the effect of lesions on slicewise connectivity values, we used linear mixed effects models in which the region of interest pair connectivity z-scores were the dependent variables. A random effect of subject was included (since connectivity values were generally higher in some subjects compared to others), and then fixed effects of subject group (control versus patient), age, gender, median slicewise temporal signal-to-noise ratio (TSNR) in grey matter, and the four binary lesion location variables (controls with all zeros). Additional analyses (not shown) including handedness as a covariate showed no significant main effects of handedness and did not appreciably affect the results presented herein. Results Connectivity in controls Connectivity results for the six region of interest pairings (where V = ventral horn, D = dorsal horn, R = right, L = left) across all slices in healthy controls (672 total slices, grey dots) are shown in Fig. 2A. Ventral-ventral (V-V) connectivity was the highest of all pairings, followed by dorsal-dorsal (D-D), and then the ipsilateral pairings. The diagonal ventral-dorsal pairings demonstrated the lowest connectivity overall. Connectivity across vertebral level (V-V and D-D) is plotted in Fig. 2B and C. There were no significant differences in V-V connectivity across vertebral levels. D-D connectivity for slices within C5 were significantly higher than the other levels obtained in our slice stack using two-tailed Wilcoxon rank-sum tests. In Fig. 2D, Z-scores for the V-V pairing in controls are grouped by gender and plotted for each subject, ordered according to subject median values. The dotted lines indicate the median z-score across all slices for each gender. Females had significantly higher V-V functional connectivity compared to males (for more on this gender effect, see the Supplementary material). Additionally, Fig. 2D demonstrates the subject-level dependence of connectivity values, motivating the linear mixed effects modelling undertaken in subsequent analyses where a random effect of subject was taken into account. Figure 2 View largeDownload slide Functional connectivity in controls and multiple sclerosis patients. (A) Box/dot plots of z-scores (i.e. connectivity strength) from every slice in controls (56 subjects × 12 slices = 672 total) for each region of interest pairing. Z-score = tanh−1(r) (dof − 3)1/2, where r is the Pearson correlation coefficient between mean timeseries from each region of interest and dof is the estimated degrees of freedom after correction for first-order autocorrelation. The mean z-score across all slices is indicated by a white midline, with black indicating the standard error of the mean and boxes indicating the standard deviation. (B and C) Box/dot plots showing z-scores for V-V and D-D across vertebral levels. Boxplots indicators are the same as in A. Two-tailed Wilcoxon rank-sum tests were performed, *P < 0.05, **P < 0.01. (D) Boxplots of z-scores for V-V in each of the 56 control subjects (12 slices/z-scores per subject), grouped by gender and ordered according to median z-score (lowest to highest within the group). The median z-score for each subject is indicated by a dark midline and boxes indicate the 25–75% quartile ranges. Whiskers extend to the most extreme z-scores not considered outliers and outliers are indicated with a black dot. Red and orange dotted lines indicate the median z-score across all males and females, respectively. (E) Group average power, after bandpass filtering. Each subject’s mean power estimated across all grey matter voxels. Wilcoxon rank-sum tests detected no significant differences between groups after division of the spectrum into 10 frequency bins. (F) Boxplots of mean z-scores for each region of interest pairing in controls (n = 56) and patients (n = 22). Subjects contributed a single mean z-score across 12 slices for each region of interest pair. Boxplots indicators are the same as previously described for D. Two-tailed Wilcoxon rank-sum tests detected no significant differences between control and patients; however, there was a trend towards decreased D-D as well as RV-LD connectivity in patients (P = 0.097 and P = 0.095, respectively). Figure 2 View largeDownload slide Functional connectivity in controls and multiple sclerosis patients. (A) Box/dot plots of z-scores (i.e. connectivity strength) from every slice in controls (56 subjects × 12 slices = 672 total) for each region of interest pairing. Z-score = tanh−1(r) (dof − 3)1/2, where r is the Pearson correlation coefficient between mean timeseries from each region of interest and dof is the estimated degrees of freedom after correction for first-order autocorrelation. The mean z-score across all slices is indicated by a white midline, with black indicating the standard error of the mean and boxes indicating the standard deviation. (B and C) Box/dot plots showing z-scores for V-V and D-D across vertebral levels. Boxplots indicators are the same as in A. Two-tailed Wilcoxon rank-sum tests were performed, *P < 0.05, **P < 0.01. (D) Boxplots of z-scores for V-V in each of the 56 control subjects (12 slices/z-scores per subject), grouped by gender and ordered according to median z-score (lowest to highest within the group). The median z-score for each subject is indicated by a dark midline and boxes indicate the 25–75% quartile ranges. Whiskers extend to the most extreme z-scores not considered outliers and outliers are indicated with a black dot. Red and orange dotted lines indicate the median z-score across all males and females, respectively. (E) Group average power, after bandpass filtering. Each subject’s mean power estimated across all grey matter voxels. Wilcoxon rank-sum tests detected no significant differences between groups after division of the spectrum into 10 frequency bins. (F) Boxplots of mean z-scores for each region of interest pairing in controls (n = 56) and patients (n = 22). Subjects contributed a single mean z-score across 12 slices for each region of interest pair. Boxplots indicators are the same as previously described for D. Two-tailed Wilcoxon rank-sum tests detected no significant differences between control and patients; however, there was a trend towards decreased D-D as well as RV-LD connectivity in patients (P = 0.097 and P = 0.095, respectively). Power spectra in patients and controls In Fig. 2E, mean power spectra for all grey matter voxels were estimated for each subject and averaged across the group for both controls (red/orange) and patients (blue), as well as an age- and gender-matched subset of the control group (12 females/10 males, mean age = 38.2 ± 10.2 years) (cyan). For this analysis, a 251-point fast Fourier transform was performed for each region of interest timecourse, yielding 100 points within the frequency range of interest (i.e. range used in bandpass filtering, 0.01–0.13 Hz). To test for differences within specific frequency bands, the range of interest was divided into 10 frequency bins with 10 points per subject. The values within each bin were averaged for each subject and Wilcoxon rank-sum tests were performed to look for group differences. No frequency bin was significantly different between the patients and either of the control groups. The qualitative increase in spectral power observed in the patient group in Fig. 2E is thus likely reflective of outliers rather than a true group difference. The spike in power observed in all groups at approximately 0.075 Hz is characteristic of a low frequency oscillation observed in brain resting state studies (Kiviniemi et al., 2005; Tong et al., 2011). It has been suggested that these spontaneous fluctuations in BOLD signal are reflective of vasomotor oscillations in arterial blood pressure/flow, e.g. Mayer-waves or C-waves (Katura et al., 2006). Fluctuations around the observed 0.075 Hz peak survived in grey matter after white matter denoising (Supplementary Fig. 1) suggesting this frequency range is related to grey matter neuronal activity. However, our white matter denoising approach does not constitute a complete suppression of signals from white matter (and surrounding spinal cord vasculature), and thus we cannot rule out the possibility that the peak at 0.075 Hz represents residual physiological noise, such as from vasomotor oscillations. There is ongoing debate about the neuronal origin of resting state BOLD fluctuations in the brain, with a recent paper suggesting that spontaneous vasomotion in the range of ∼0.1 Hz tracks the envelope of gamma-band electrical activity in the cortex (Mateo et al., 2017), while another suggests that resting state BOLD fluctuations are influenced by unrelated haemodynamic processes and only weakly correlate with neural activity (Winder et al., 2017). Future research is required to further investigate the source of spinal cord BOLD signals and their potential neuronal origins, such as via intraspinal electrical recording. Further analyses assessing the contribution of narrower frequency bands to our functional connectivity metrics are provided in Supplementary Fig. 2. Figure 3 View largeDownload slide Lesion ratings. Patient lesion ratings were conducted by three independent raters and involved determination for each slice whether a lesion was present in ventral, dorsal, left lateral, and/or right lateral column of the spinal cord white matter. Final consensus ratings were determined and are depicted for each subject. Example slices, which were rated as containing a lesion, are provided. The straightened spinal cord and representative cervical slice images to the right were adapted from recently developed spinal cord templates included in the Spinal Cord Toolbox software package (Fonov et al., 2014; De Leener et al., 2017, 2018). Figure 3 View largeDownload slide Lesion ratings. Patient lesion ratings were conducted by three independent raters and involved determination for each slice whether a lesion was present in ventral, dorsal, left lateral, and/or right lateral column of the spinal cord white matter. Final consensus ratings were determined and are depicted for each subject. Example slices, which were rated as containing a lesion, are provided. The straightened spinal cord and representative cervical slice images to the right were adapted from recently developed spinal cord templates included in the Spinal Cord Toolbox software package (Fonov et al., 2014; De Leener et al., 2017, 2018). Comparison of functional connectivity in patients to controls To look for differences in functional connectivity between patient and control groups, the mean z-score across 12 slices for each subject/region of interest pairing was first computed (Fig. 2F). Tests of difference were conducted on the mean values using two-tailed Wilcoxon rank-sum tests for each region of interest pairing. There were no significant differences in z-scores between subject groups for any of the region of interest pairs, although a trend towards decreased D-D as well as RV-LD connectivity was observed in patients (P = 0.097 and P = 0.095, respectively). While the rank-sum tests served as a first pass for looking at group differences in the raw connectivity values, they did not control for other variables such as the significant difference in age between the patient and control groups or potential gender effects, as is demonstrated in Fig. 2D. To further assess group differences while controlling for other variables, linear models were fit using multiple regression. Four predictor variables were included in each model, group (Control/Patient), age, gender (M/F), and TSNR (calculated after white matter denoising but prior to bandpass filtering). The main effect of each variable in predicting the average region of interest pair z-score over 12 slices is depicted in Supplementary Fig. 4 along with the confidence intervals for these values. None of the region of interest pairings demonstrated a significant main effect of group. One challenge in interpreting these results is that our patient sample was significantly older, potentially conflating effects of disease and normal ageing, for instance in the D-D connectivity. Finally, we assessed correlations between clinical measures of disability and average connectivity for each region of interest pairing in the patient group. We found that no disability measure [i.e. Expanded Disability Status Scale (EDSS), disease duration, and timed 25-foot walk test] was significantly directly correlated with mean connectivity strength in patients (Supplementary Fig. 3). When controlling for nuisance variables (i.e. age, gender, and TSNR), however, we found modest positive associations of EDSS and 25-foot walk times with D-D connectivity, as well a positive association between 25-foot walk times and V-V connectivity (Supplementary Fig. 5). These findings indicate that increasing disability is related to increased bilateral connectivity, and it is possible that this reflects a compensatory mechanism among bilateral spinal cord networks. However, these observations should be considered preliminary and more work with a larger sample size is necessary to confirm this finding. Lesion ratings analysis To assess the level of agreement between the three independent raters, we looked at the per cent agreement across all possible lesions locations (i.e. lesion present or not) in which both raters provided a rating (note that raters had the option to not rate an entire slice due to artefact). Between each pair of raters, the per cent agreement was 90.7% (1028 total locations), 88.2% (952 total locations), and 91.4% (940 total locations), respectively. Following this analysis, a consensus meeting was held in which the discrepancies were reviewed and final ratings were agreed upon. In five slices (out of 264, ∼2%), a rating was not assigned because of artefacts. Of the possible 1036 locations considered for the final analysis, 265 locations were rated as containing a lesion. Overall, 19 of 22 patients presented with a lesion somewhere within the imaged volume. The final consensus lesion ratings were taken as the ground truth for subsequent analyses and are visually depicted in Fig. 3. Local effects of lesions on functional connectivity To assess the effect of lesions on connectivity within a slice (local), linear mixed effects models were fit individually for each of the region of interest pairings, with the slicewise connectivity z-score as the outcome variable, a random effect of subject, and the remaining variables as fixed effects. Control data were included in these models, with lesion values at zero. The relevant statistical results for the six models are provided in Table 2, and the full statistical data including parameter estimates for nuisance variables are provided in Supplementary Table 2. Effect sizes and polarity of the estimate (whether positive or negative) for each lesion location are visually depicted in Fig. 4. Ventral column lesions demonstrated a significant negative association with local V-V connectivity such that the presence of a ventral lesion was associated with decreased V-V connectivity after controlling for all other effects, including the effects of lesions in other locations. In the D-D model, the presence of a left lateral lesion demonstrated a significant positive association with D-D connectivity. Additionally, higher RV-RD connectivity was associated with the presence of right lateral lesions and higher LV-RD connectivity was associated with dorsal column lesions. Table 2 Local effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Table 2 Local effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Figure 4 View largeDownload slide Local effects of lesions on functional connectivity. Visual representation of effect sizes for lesion presence (ventral, dorsal, left lateral, and right lateral column) in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 2. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue circles indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the column, which they represent. All patient (n = 22) and control (n = 56) data were included in the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for lesion values. Full statistical results are provided in Supplementary Table 2. Note that all lesion variables were included in the same model, such that the estimates provided are the independent effect of a lesion within the column on connectivity after accounting for all other effects, including lesions in other columns. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes in table indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. FC = functional connectivity. Figure 4 View largeDownload slide Local effects of lesions on functional connectivity. Visual representation of effect sizes for lesion presence (ventral, dorsal, left lateral, and right lateral column) in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 2. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue circles indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the column, which they represent. All patient (n = 22) and control (n = 56) data were included in the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for lesion values. Full statistical results are provided in Supplementary Table 2. Note that all lesion variables were included in the same model, such that the estimates provided are the independent effect of a lesion within the column on connectivity after accounting for all other effects, including lesions in other columns. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes in table indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. FC = functional connectivity. Distal effects of lesions on functional connectivity To assess the effect of lesions on connectivity above and below the lesion, we first classified lesions as occurring in ascending and/or descending white matter tracts (see schematic in Fig. 5). Dorsal column lesions were classified as ascending tract lesions. Ventral column lesions were classified as descending tract lesions. Since both ascending and descending tracts are found in the lateral columns, these lesions were counted for both classes. Our anatomical imaging was not of sufficient resolution to confidently distinguish lesions within specific lateral column tracts. Furthermore, many of the lateral lesions encompassed large areas of the column, thus likely affecting both ascending and descending tracts (see lesion examples in Fig. 3). Slices were then assigned a value of 1 if they occurred either above an ascending tract lesion or below a descending tract lesion. Linear mixed effects models were then individually fit for each of the region of interest pairs and the relevant results are provided in Table 3. The full statistical results including parameter estimates for nuisance variables are provided in Supplementary Table 3. A visual depiction of the effect sizes for lesion location are displayed in Fig. 5. Though not significant at P < 0.05, a modest negative effect was observed for V-V connectivity in slices below a descending tract lesion after controlling for other variables (P = 0.131). LD-RV connectivity in slices below a descending tract lesion was significantly decreased (P = 0.011). Within the same model, an increase in LD-RV connectivity in slices above an ascending tract lesion approached significance (P = 0.057). Table 3 Distal effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Table 3 Distal effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Figure 5 View largeDownload slide Distal effects of lesions on functional connectivity. Visual representation of effect sizes for distal presence of lesion in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 3. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue arrows indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the regions which they represent. Slices were rated as containing a descending tract lesion if a lesion was present in at least one of the ventral, left lateral, or right lateral tracts; and a descending tract lesion if a lesion was present in at least one of the dorsal, left lateral, or right lateral tracts. Based on this classification, two independent variables were constructed in which slices were assigned a binary value indicating whether they did (1) or did not (0) occur (i) above (superior to) an ascending tract lesion (‘Above’ variable in Table 3, upward facing arrows in figure); and (ii) below (inferior to) a descending tract lesion (‘Below’ variable in Table, downward facing arrows in figure). Slices could receive only binary classifications of 1 or 0 for the Above and Below variables and thus were not additive, for instance in cases where a slice was above multiple ascending tract lesions. All patient (n = 22) and control (n = 56) data were included the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for Above and Below values. Full statistical results are provided in Supplementary Table 3. Note that the Above and Below variables were included in the same model, such that the estimates provided are the independent effect of a slice occurring above an ascending tract lesion after accounting for all other effects, including the slice occurring below a descending tract lesion. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. Figure 5 View largeDownload slide Distal effects of lesions on functional connectivity. Visual representation of effect sizes for distal presence of lesion in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 3. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue arrows indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the regions which they represent. Slices were rated as containing a descending tract lesion if a lesion was present in at least one of the ventral, left lateral, or right lateral tracts; and a descending tract lesion if a lesion was present in at least one of the dorsal, left lateral, or right lateral tracts. Based on this classification, two independent variables were constructed in which slices were assigned a binary value indicating whether they did (1) or did not (0) occur (i) above (superior to) an ascending tract lesion (‘Above’ variable in Table 3, upward facing arrows in figure); and (ii) below (inferior to) a descending tract lesion (‘Below’ variable in Table, downward facing arrows in figure). Slices could receive only binary classifications of 1 or 0 for the Above and Below variables and thus were not additive, for instance in cases where a slice was above multiple ascending tract lesions. All patient (n = 22) and control (n = 56) data were included the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for Above and Below values. Full statistical results are provided in Supplementary Table 3. Note that the Above and Below variables were included in the same model, such that the estimates provided are the independent effect of a slice occurring above an ascending tract lesion after accounting for all other effects, including the slice occurring below a descending tract lesion. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. Discussion The present study provides the first assessment of functional MRI-derived functional connectivity in the cord of a diseased population. This builds upon recent reports of intrinsic resting state BOLD activity and correlation in the healthy human spinal cord. Spinal cord functional MRI presents unique challenges, such as requiring higher resolution than studies in the brain due to the small structures of interest. We took advantage of ultra-high field 7 T imaging, which provides increased BOLD contrast as well as increased sensitivity in the detection of multiple sclerosis lesions. Our previous work provided a methodological framework for the measurement of functional connectivity in patients, including a 3D-multishot functional MRI sequence, which reduced susceptibility distortions and an optimized preprocessing and denoising pipeline. The sample of healthy controls considered in this analysis represents the largest described in the spinal cord functional MRI literature to date and further supports the existence of distinct and robust functional networks within human cervical cord grey matter. In both patients and controls, the strongest correlations were observed among bilateral ventral grey matter regions, followed by those among dorsal grey matter regions. Importantly, while more research is needed to confirm the neuronal origin of spinal cord BOLD signals, our work suggests spinal cord networks are detectable in the absence of a particular task or stimulation, making the acquisition and translation of these measures more feasible for future clinical application. We found no significant differences in mean connectivity between multiple sclerosis patients and controls at the group level (Fig. 2F), but did detect both decreased and increased connectivity among grey matter regions in relation to local multiple sclerosis lesions (Fig. 4 and Table 2). For instance, we found a significant decrease in V-V connectivity in the presence of ventral column lesions, suggesting that damage to the ventral column white matter may be particularly detrimental to motor network synchronization within the spinal cord. It is possible that ventral column lesions characterize a direct disconnection between intra-segmental ventral grey matter circuits. For other region pairs in which we found increased connectivity near lesion sites (D-D with L lesion, RV-RD with R lesion and trend with D lesion, and LV-RD with D lesion) one hypothesis is that increased connectivity in the presence of lesions and/or axonal loss (Kearney et al., 2015) reflects a compensatory mechanism within local grey matter circuits. One study by Agosta et al. (2008b) found that there was increased bilateral recruitment of grey matter regions in patients with multiple sclerosis in response to stimulation, and that an increased signal change to proprioceptive stimulation was associated with reduced fractional anisotropy in the cord, suggesting that increasing white matter damage is related to a compensatory increase in activity and/or metabolic demand in functioning grey matter. Furthermore, the present findings are in-line with converging evidence of compensatory increases in brain connectivity in relapsing-remitting multiple sclerosis in sensorimotor network as well as other resting state networks, particularly those showing increased connectivity with increased structural damage. For instance, a study using resting functional MRI and diffusion tensor imaging in minimally disabled relapsing-remitting patients compared to controls found that regions of the default mode network were functionally overconnected while structurally disconnected, due to damage in long-range white matter tracts, suggesting a compensatory mechanism of a network, albeit at a more global scale (Zhou et al., 2014). Our findings of increased local connectivity for some region pairs/lesion sites may too reflect a compensatory effect of white matter damage. Though not mutually exclusive, an alternative mechanism has been proposed to explain overactivation and increased bilateral recruitment of cervical cord grey matter in multiple sclerosis, which involves a disruption of inhibitory spinal interneurons, such that there is failure to inhibit a widespread network of primary spinal neurons. Histopathological assessment of multiple sclerosis spinal cord tissue has demonstrated both a reduced number and decreased cell size of interneurons in the upper cervical levels compared to control tissue (Gilmore et al., 2009). An abnormal increase in the intrinsic coactivation (i.e. connectivity) of grey matter in a resting state may also be reflective of interneuron dysfunction. If so, our findings may indicate interneuron effects on intrinsic activity are most pronounced near lesion sites. While further work is necessary to distinguish between these hypotheses, our findings suggest that the functional consequences of multiple sclerosis cord lesions as measured through fluctuations in grey matter BOLD signals, such as from compensatory mechanisms and/or interneuron abnormalities, are detectable even at rest. Given the connectivity strength was relatively low between ipsilateral (LV-LD and RV-RD) and diagonal (LV-RD and RV-LD) region of interest pairs as compared to V-V and D-D (Fig. 2A), caution should be taken in interpreting significant effects involving these pairs (i.e. Table 2, RV-RD/R lesion, LV-RD/D lesion, and Table 3, LD-RV/Below). It may be, however, that synchronous activity between ipsilateral and diagonal grey matter horns occurs in relatively transient intervals at rest, such that the static connectivity over the entire imaging run is low overall. If so, significant effects involving dorsal-ventral horn connectivity may be indicative of relevant physiological processes. Future work using dynamic connectivity analyses, for instance, may shed light on the degree of transient synchrony between ipsilateral and diagonal regions of grey matter in the cord. In terms of the distal effects of lesions on connectivity, a trend was observed in which V-V connectivity was decreased when slices were below a lesion within a descending tract (Fig. 5, P = 0.131). Given ventral grey matter is largely composed of motor neurons that receive input from descending corticospinal axons, it seems reasonable that lesions within these tracts disrupt synchronization of motor neuron populations below the site of pathology. Furthermore, our findings of local and distal effects of structural damage on V-V connectivity are in agreement with the spinal cord injury study by Chen et al. (2015), which showed decreased connectivity within and below the site of injury, but not above. It is important to note, however, that the pathophysiological course of acute spinal cord injury and its subsequent effects on functional connectivity are likely to be qualitatively different from that of demyelinating lesions in multiple sclerosis. The cellular and biochemical cascades accompanied by, for example, gross inflammation and/or breakdown of blood–spinal cord barrier in response to acute injury may have more generalized effects on cord function inferior to injury whereas multiple sclerosis lesion effects, especially during periods of remission, may be more tract-specific (Hausmann, 2003). Nevertheless, the prospects for non-invasive measurement of cord functional connectivity in spinal cord injury and multiple sclerosis patient populations provide an interesting avenue of research for understanding both similarities and differences among dysfunctional cord pathologies. There was a general lack of group level differences in mean functional connectivity metrics between patients with multiple sclerosis and healthy controls (Fig. 2F), especially after controlling for confounding factors (Supplementary Fig. 4). On one hand, it is perhaps encouraging that our connectivity measures are stable across groups and region of interest pairs. However, there are several reasons we may have failed to detect a group difference in mean connectivity values. Our patient sample was relatively heterogeneous in terms of lesion load, clinical disability, and disease duration. This could have potentially conflated early disease stage compensatory increases in connectivity with later stage decreases across the patient group, such that any differences compared to controls were washed out. It is also possible that disease effects could not be teased apart from age effects in our statistical models, given our patient sample was older than controls. Another possibility is that intrinsic spinal cord activity is affected at more acute time windows in relapsing-remitting multiple sclerosis, such as during a relapse phase, and recovers to more normal levels in remitting phases. This would be in-line with both Droby et al. (2016) and Chen et al. (2015), who observed normalization of functional connectivity within several months after acute lesion appearance in the brain and induced injury in cord, respectively. In contrast, the observations from functional MRI activation studies in multiple sclerosis seem somewhat counter to the present findings, where group-level overactivation effects were observed even in relapsing-remitting cohorts. It is important to note that activation in response to a stimulus is not necessarily related to functional connectivity strength, which is instead a measure of the relative synchrony of activity between two regions over time. While functional connectivity patterns at rest are largely analogous to those observed in overt task-based activation studies, resting state networks in the brain are found to be more bilateral in nature. Our results demonstrating distinct, bilateral connectivity between ventral grey matter and between dorsal grey matter in the cord corroborate the general finding of symmetrical networks in the brain. A final caveat to be considered in interpretation of this work as compared to the aforementioned studies of cord activation in multiple sclerosis is that all of these previous studies were performed on 1.5 T scanners and used spin echo functional MRI acquisitions, which are much less sensitive to traditional T2*-weighted BOLD effects and instead are suggested to measure a related, but alternative phenomenon termed SEEP (signal enhancement by extravascular protons) (Stroman et al., 2002, 2003). While the benefits of SEEP functional MRI in the spinal cord have been described (Stroman et al., 2005, 2014), this contrast mechanism still remains controversial among the wider functional MRI community in its ability to measure effects specific to neuronal activation (Jochimsen et al., 2005; Bouwman et al., 2008; Summers and Brooks, 2014). More recent work looking at stimulus-evoked activity and resting state fluctuations in cord, including the present study, have largely transitioned to the employment of BOLD-sensitive acquisitions. Going forward, it will be important to confirm the findings of multiple sclerosis-related overactivation in response to stimulation using BOLD functional MRI. Limitations The present study provides a novel assessment of spinal cord functional connectivity in patients with multiple sclerosis compared to healthy controls, including observations of differential effects of column-specific lesion location on synchronous intrinsic BOLD fluctuations in cervical cord grey matter, and there are a number of limitations to consider in interpreting this work (for a more thorough discussion of some of these limitations, please see the Supplementary material): (i) our rating strategy only provided a coarse identification of lesion location, at the level of columns. Damage to tracts, which are known to specifically innervate cervical cord grey matter (e.g. lateral corticospinal tract), may have more robust effects on connectivity of these regions; (ii) while we took a region of interest-based approach dividing spinal cord grey matter into four horns, connectivity in the cord may be divisible into smaller subnetworks and could benefit from analyses using smaller regions of interest or spatial independent component analysis, which looks for patterns at the voxel level (Kong et al., 2014); (iii) grey matter lesions exist in the spinal cord in multiple sclerosis (Gilmore et al., 2009), as has been appreciated in the brain (Calabrese et al., 2013; Harrison et al., 2015). Future studies should adopt imaging strategies to distinguish healthy and lesioned grey matter in the spinal cord, as this damage may influence connectivity measures (Nelson et al., 2008; De Graaf et al., 2012; Sinnecker et al., 2012); (iv) an important limitation is that, while we looked at the effects of relatively local cervical lesions, structural damage located outside of our imaged volume could have influenced connectivity of the cervical regions analysed in this study. Brain lesions affecting corticospinal tracts, for instance, may modulate spontaneous ventral horn activity or connectivity. With the cervical spine RF-coil used in this study, we were not able to acquire whole-brain images to assess brain lesion burden or location. Future work is needed to understand how brain lesions influence spinal cord function in multiple sclerosis, and a more comprehensive assessment of CNS abnormalities (e.g. inclusion of whole-brain imaging) will be necessary to distinguish upstream/downstream contributions to spinal cord connectivity patterns; (v) it is possible that muscle relaxants and/or sedative drugs such as baclofen or benzodiazepine, which are commonly used to treat patients with spinal cord disorders, have an effect on spinal cord network connectivity. In our patient sample, 7 of 22 subjects were prescribed one of these medications (Supplementary Table 1). We assessed the effect of these medications on mean connectivity and found no significant associations with any region of interest pairings (Supplementary Fig. 5). However, future work with a larger sample size is necessary to confirm this observation and researchers may consider use of these medications as exclusion criteria; and (vi) a final limitation concerns the current practicality of acquiring spinal functional connectivity measurements in a clinical setting, given the relatively complex image processing pipeline used here, which requires some expertise as well as use of a 7 T MRI machine, which are still not widely available. The spinal cord MRI community continues to grow, and along with it is continued development of more user-friendly processing software, similar to what is available for brain imaging. The Spinal Cord Toolbox provides the first comprehensive software package for the processing of spinal cord MRI data, including functional MRI analyses (De Leener et al., 2017). Strategies for robust automatic segmentation of the spinal cord have recently been described as well, which reduce processing time and provide further promise for future clinical application (Asman et al., 2014; De Leener et al., 2016; Prados et al., 2016, 2017; Xu et al., 2016). Regarding 7 T MRI, 510(k) approval for the Siemens 7 T Magnetom Terra system was recently announced by the US Food and Drug Administration, and it is anticipated that this announcement will, in time, facilitate more widespread adoption of 7 T MRI. Nevertheless, several groups have now shown spinal cord functional connectivity is observable at 3 T (e.g. Eippert et al., 2017b; Harita and Stroman, 2017; Weber et al., 2018), suggesting these measures can feasibly be acquired on current clinical scanners. Acknowledgement The authors thank Avyay Reddy for performing quality assurance on preprocessed data. Funding This research was supported by the following funding sources: Novartis IIRP-1456/VUMC 41030 (Pawate), NIH/NIBIB 4T32EB014841-04 (Gore), NIH/NINDS RO1NS092961 (Gore), NIH/NINDS R01NS104149 (Gore), DOD W81XWH-13-1-0073 (Smith), National MS Society, and NIH/NIBIB K99/R00 EB016689 (Barry). The project was also supported by the National Center for Research Resources, grant UL1RR024975-01, which is now at the National Center for Advancing Translational Sciences, grant 2UL1TR000445-06. Supplementary material Supplementary material is available at Brain online. Abbreviation Abbreviation BOLD blood oxygenation level-dependent References Agosta F, Valsasina P, Absinta M, Sala S, Caputo D, Filippi M. Primary progressive multiple sclerosis: tactile-associated functional MR activity in the cervical spinal cord. Radiology  2009; 253: 209– 15. Google Scholar CrossRef Search ADS PubMed  Agosta F, Valsasina P, Caputo D, Stroman PW, Filippi M. Tactile-associated recruitment of the cervical cord is altered in patients with multiple sclerosis. Neuroimage  2008a; 39: 1542– 8. Google Scholar CrossRef Search ADS   Agosta F, Valsasina P, Rocca MA, Caputo D, Sala S, Judica E, et al.   Evidence for enhanced functional activity of cervical cord in relapsing multiple sclerosis. Magn Reson Med  2008b; 59: 1035– 42. Google Scholar CrossRef Search ADS   Asman AJ, Bryan FW, Smith SA, Reich DS, Landman BA. Groupwise multi-atlas segmentation of the spinal cord’s internal structure. Med Image Anal  2014; 18: 460– 71. Google Scholar CrossRef Search ADS PubMed  Barry RL, Gore JC, Smith SA. Development of BOLD-sensitive protocols for imaging the human spinal cord at 7 Tesla. Proc Intl Soc Mag Reson Med  2013; 21: 2278. Barry RL, Rogers BP, Conrad BN, Smith SA, Gore JC. Reproducibility of resting state spinal cord networks in healthy volunteers at 7 Tesla. Neuroimage  2016; 133: 31– 40. Google Scholar CrossRef Search ADS PubMed  Barry RL, Smith SA, Dula AN, Gore JC. Resting state functional connectivity in the human spinal cord. Elife  2014; 3: e02812. Google Scholar PubMed  Barry RL, Strother SC, Gatenby JC, Gore JC. Data-driven optimization and evaluation of 2D EPI and 3D PRESTO for BOLD fMRI at 7 Tesla: I. Focal coverage. Neuroimage  2011; 55: 1034– 43. Google Scholar CrossRef Search ADS PubMed  Barry RL, Vannesjo SJ, By S, Gore JC, Smith SA. Spinal cord MRI at 7T. Neuroimage  2018; 168: 437– 51. Google Scholar CrossRef Search ADS PubMed  Bodini B, Battaglini M, De Stefano N, Khaleeli Z, Barkhof F, Chard D, et al.   T2 lesion location really matters: a 10 year follow-up study in primary progressive multiple sclerosis. J Neurol Neurosurg Psychiatry  2011; 82: 72– 7. Google Scholar CrossRef Search ADS PubMed  Bouwman CJC, Wilmink JT, Mess WH, Backes WH. Spinal cord functional MRI at 3 T: gradient echo echo-planar imaging versus turbo spin echo. Neuroimage  2008; 43: 288– 96. Google Scholar CrossRef Search ADS PubMed  Brown TG. The intrinsic factors in the act of progression in the mammal. Proc R Soc London Ser B  1911; 84: 308– 19. Google Scholar CrossRef Search ADS   Cadotte DW, Stroman PW, Mikulis D, Fehlings MG. A systematic review of spinal fMRI research: outlining the elements of experimental design. J Neurosurgery Spine  2012; 17: 102– 18. Google Scholar CrossRef Search ADS   Calabrese M, Favaretto A, Martini V, Gallo P. Grey matter lesions in MS: from histology to clinical implications. Prion  2013; 7: 20– 7. Google Scholar CrossRef Search ADS PubMed  Chen LM, Mishra A, Yang PF, Wang F, Gore JC. Injury alters intrinsic functional connectivity within the primate spinal cord. Proc Natl Acad Sci USA  2015; 112: 5991– 6. Google Scholar CrossRef Search ADS PubMed  Cohen-Adad J, Wheeler-Kingshott C. Quantitative MRI of the Spinal Cord . New York: Academic Press; 2014. Conrad BN, Maki S, Watchmaker JM, Bailey BA, Barry RL, Smith SA, et al.   BOLD Signal changes in spinal cord with hypercapnia. Proc Intl Soc Mag Reson Med  2017; 25: 5319. Cox RW. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res  1996; 29: 162– 73. Google Scholar CrossRef Search ADS PubMed  De Leener B, Fonov VS, Collins DL, Callot V, Stikov N, Cohen-Adad J. PAM50: unbiased multimodal template of the brainstem and spinal cord aligned with the ICBM152 space. Neuroimage  2018; 165: 170– 9. Google Scholar CrossRef Search ADS PubMed  De Leener B, Lévy S, Dupont SM, Fonov VS, Stikov N, Louis Collins D, et al.   SCT: Spinal Cord Toolbox, an open-source software for processing spinal cord MRI data. Neuroimage  2017; 145: 24– 43. Google Scholar CrossRef Search ADS PubMed  De Leener B, Taso M, Cohen-Adad J, Callot V. Segmentation of the human spinal cord. MAGMA  2016; 29: 125– 53. Google Scholar CrossRef Search ADS PubMed  Droby A, Yuen KSL, Muthuraman M, Reitz SC, Fleischer V, Klein J, et al.   Changes in brain functional connectivity patterns are driven by an individual lesion in MS: a resting-state fMRI study. Brain Imaging Behav  2016; 10: 1117– 26. Google Scholar CrossRef Search ADS PubMed  Dula AN, Pawate S, Dortch RD, Barry RL, George-Durrett KM, Lyttle BD, et al.   Magnetic resonance imaging of the cervical spinal cord in multiple sclerosis at 7T. Mult Scler  2016; 22: 320– 8. Google Scholar CrossRef Search ADS PubMed  Eippert F, Kong Y, Jenkinson M, Tracey I, Brooks JCW. Denoising spinal cord fMRI data: approaches to acquisition and analysis. Neuroimage  2017a; 154: 255– 66. Google Scholar CrossRef Search ADS   Eippert F, Kong Y, Winkler AM, Andersson JL, Finsterbusch J, Büchel C, et al.   Investigating resting-state functional connectivity in the cervical spinal cord at 3T. Neuroimage  2017b; 147: 589– 601. Google Scholar CrossRef Search ADS   Faivre A, Rico A, Zaaraoui W, Crespy L, Reuter F, Wybrecht D, et al.   Assessing brain connectivity at rest is clinically relevant in early multiple sclerosis. Mult Scler J  2012; 18: 1251– 8. Google Scholar CrossRef Search ADS   Filippi M, Rocca MA. Present and future of fMRI in multiple sclerosis. Expert Rev Neurother  2013; 13: 27– 31. Google Scholar CrossRef Search ADS PubMed  Fonov VS, Le Troter A, Taso M, De Leener B, Lévêque G, Benhamou M, et al.   Framework for integrated MRI average of the spinal cord white and gray matter: the MNI-Poly-AMU template. Neuroimage  2014; 102: 817– 27. Google Scholar CrossRef Search ADS PubMed  Gilmore CP, Deluca GC, Bö L, Owens T, Lowe J, Esiri MM, et al.   Spinal cord neuronal pathology in multiple sclerosis. Brain Pathol  2009; 19: 642– 9. Google Scholar CrossRef Search ADS PubMed  Giove F, Garreffa G, Giulietti G, Mangia S, Colonnese C, Maraviglia B. Issues about the fMRI of the human spinal cord. Magn Reson Imaging  2004; 22: 1505– 16. Google Scholar CrossRef Search ADS PubMed  Glover GH, Li TQ, Ress D. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Reson Med  2000; 44: 162– 7. Google Scholar CrossRef Search ADS PubMed  De Graaf WL, Zwanenburg JJM, Visser F, Wattjes MP, Pouwels PJW, Geurts JJG, et al.   Lesion detection at seven Tesla in multiple sclerosis using magnetisation prepared 3D-FLAIR and 3D-DIR. Eur Radiol  2012; 22: 221– 31. Google Scholar CrossRef Search ADS PubMed  Harita S, Stroman PW. Confirmation of resting-state BOLD fluctuations in the human brainstem and spinal cord after identification and removal of physiological noise. Magn Reson Med  2017; 78: 2149– 56. Google Scholar CrossRef Search ADS PubMed  Harrison DM, Roy S, Oh J, Izbudak I, Pham D, Courtney S, et al.   Association of cortical lesion burden on 7-T magnetic resonance imaging with cognition and disability in multiple sclerosis. JAMA Neurol  2015; 72: 1004– 12. Google Scholar CrossRef Search ADS PubMed  Hausmann ON. Post-traumatic inflammation following spinal cord injury. Spinal Cord  2003; 41: 369– 78. Google Scholar CrossRef Search ADS PubMed  Jochimsen TH, Norris DG, Möller HE. Is there a change in water proton density associated with functional magnetic resonance imaging? Magn Reson Med  2005; 53: 470– 3. Google Scholar CrossRef Search ADS PubMed  Kandel ER, Schwartz JH, Jessell TM, Siegelbaum SA, Hudspeth AJ, Mack S. Principles of neural science . 5th edn. New York: McGraw-Hill Companies, Incorporated; 2013. Katura T, Tanaka N, Obata A, Sato H, Maki A. Quantitative evaluation of interrelations between spontaneous low-frequency oscillations in cerebral hemodynamics and systemic cardiovascular dynamics. Neuroimage  2006; 31: 1592– 600. Google Scholar CrossRef Search ADS PubMed  Kearney H, Miller DH, Ciccarelli O. Spinal cord MRI in multiple sclerosis-diagnostic, prognostic and clinical value. Nat Rev Neurol  2015; 11: 327– 38. Google Scholar CrossRef Search ADS PubMed  Kiviniemi V, Ruohonen J, Tervonen O. Separation of physiological very low frequency fluctuation from aliasing by switched sampling interval fMRI scans. Magn Reson Imaging  2005; 23: 41– 6. Google Scholar CrossRef Search ADS PubMed  Kong Y, Eippert F, Beckmann CF, Andersson J, Finsterbusch J, Büchel C, et al.   Intrinsically organized resting state networks in the human spinal cord. Proc Natl Acad Sci USA  2014; 111: 18067– 72. Google Scholar CrossRef Search ADS PubMed  Liu X, Zhou F, Li X, Qian W, Cui J, Zhou IY, et al.   Organization of the intrinsic functional network in the cervical spinal cord: a resting state functional MRI study. Neuroscience  2016; 336: 30– 8. Google Scholar CrossRef Search ADS PubMed  Maieron M, Iannetti GD, Bodurka J, Tracey I, Bandettini PA, Porro CA. Functional responses in the human spinal cord during willed motor actions: evidence for side- and rate-dependent activity. J Neurosci  2007; 27: 4182– 90. Google Scholar CrossRef Search ADS PubMed  Mateo C, Knutsen PM, Tsai PS, Shih AY, Kleinfeld D. Entrainment of arteriole vasomotor fluctuations by neural activity is a basis of blood-oxygenation-level-dependent ‘resting-state’ connectivity. Neuron  2017; 96: 936– 48.e3. Google Scholar CrossRef Search ADS PubMed  McCrea DA, Rybak IA. Organization of mammalian locomotor rhythm and pattern generation. Brain Res Rev  2008; 57: 134– 46. Google Scholar CrossRef Search ADS PubMed  Nelson F, Poonawalla A, Hou P, Wolinsky JS, Narayana PA. 3D MPRAGE improves classification of cortical lesions in multiple sclerosis. Mult Scler  2008; 14: 1214– 19. Google Scholar CrossRef Search ADS PubMed  Prados F, Ashburner J, Blaiotta C, Brosch T, Carballido-Gamio J, Cardoso MJ, et al.   Spinal cord grey matter segmentation challenge. Neuroimage  2017; 152: 312– 29. Google Scholar CrossRef Search ADS PubMed  Prados F, Cardoso MJ, Yiannakas MC, Hoy LR, Tebaldi E, Kearney H, et al.   Fully automated grey and white matter spinal cord segmentation. Sci Rep  2016; 6: 36151. Google Scholar CrossRef Search ADS PubMed  Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: Sensitivity encoding for fast MRI. Magn Reson Med  1999; 42: 952– 62. Google Scholar CrossRef Search ADS PubMed  Purves D, Augustine G, Fitzpatrick D, Hall W, Lamantia AS, White L. Neuroscience . 5th edn. Sunderland, MA: Sinauer Associates, Incorporated; 2012. Rocca MA, Absinta M, Valsasina P, Copetti M, Caputo D, Comi G, et al.   Abnormal cervical cord function contributes to fatigue in multiple sclerosis. Mult Scler J  2012a; 18: 1552– 9. Google Scholar CrossRef Search ADS   Rocca MA, Amato MP, De Stefano N, Enzinger C, Geurts JJ, Penner IK, et al.   Clinical and imaging assessment of cognitive dysfunction in multiple sclerosis. Lancet Neurol  2015; 14: 302– 417. Google Scholar CrossRef Search ADS PubMed  Rocca MA, Valsasina P, Martinelli V, Misci P, Falini A, Comi G, et al.   Large-scale neuronal network dysfunction in relapsing-remitting multiple sclerosis. Neurology  2012b; 79: 1449– 57. Google Scholar CrossRef Search ADS   Rogers BP, Gore JC. Empirical comparison of sources of variation for FMRI connectivity analysis. PLoS One  2008; 3: e3708. Google Scholar CrossRef Search ADS PubMed  Roosendaal SD, Schoonheim MM, Hulst HE, Sanz-Arigita EJ, Smith SM, Geurts JJG, et al.   Resting state networks change in clinically isolated syndrome. Brain  2010; 133: 1612– 21. Google Scholar CrossRef Search ADS PubMed  Sacco R, Bonavita S, Esposito F, Tedeschi G, Gallo A. The contribution of resting state networks to the study of cortical reorganization in MS. Mult Scler Int  2013; 2013: 857807. Google Scholar PubMed  San Emeterio Nateras O, Yu F, Muir ER, Bazan C, Franklin CG, Li W, et al.   Intrinsic resting-state functional connectivity in the human spinal cord at 3.0 T. Radiology  2016; 279: 262– 8. Google Scholar CrossRef Search ADS PubMed  Sinnecker T, Dörr J, Pfueller CF, Harms L, Ruprecht K, Jarius S, et al.   Distinct lesion morphology at 7-T MRI differentiates neuromyelitis optica from multiple sclerosis. Neurology  2012; 79: 708– 14. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Bosma RL, Tsyben A. Somatotopic arrangement of thermal sensory regions in the healthy human spinal cord determined by means of spinal cord functional MRI. Magn Reson Med  2012; 68: 923– 31. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Kornelsen J, Lawrence J, Malisza KL. Functional magnetic resonance imaging based on SEEP contrast: response function and anatomical specificity. Magn Reson Imaging  2005; 23: 843– 50. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Krause V, Malisza KL, Frankenstein UN, Tomanek B. Extravascular proton-density changes as a non-BOLD component of contrast in fMRI of the human spinal cord. Magn Reson Med  2002; 48: 122– 7. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Tomanek B, Krause V, Frankenstein UN, Malisza KL. Functional magnetic resonance imaging of the human brain based on signal enhancement by extravascular protons (SEEP fMRI). Magn Reson Med  2003; 49: 433– 9. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Wheeler-Kingshott C, Bacon M, Schwab JM, Bosma R, Brooks J, et al.   The current state-of-the-art of spinal cord imaging: methods. Neuroimage  2014; 84: 1070– 81. Google Scholar CrossRef Search ADS PubMed  Summers PE, Brooks JCW. Chapter 4.1—Spinal cord fMRI. In: Quantitative MRI of the spinal cord . New York: Academic Press; 2014. p. 221– 39. Tong Y, Lindsey KP, deB Frederick B. Partitioning of physiological noise signals in the brain with concurrent near-infrared spectroscopy and fMRI. J Cereb Blood Flow Metab  2011; 31: 2352– 62. Google Scholar CrossRef Search ADS PubMed  Valsasina P, Agosta F, Absinta M, Sala S, Caputo D, Filippi M. Cervical cord functional MRI changes in relapse-onset MS patients. J Neurol Neurosurg Psychiatry  2010; 81: 405– 8. Google Scholar CrossRef Search ADS PubMed  Valsasina P, Rocca MA, Absinta M, Agosta F, Caputo D, Comi G, et al.   Cervical cord FMRI abnormalities differ between the progressive forms of multiple sclerosis. Hum Brain Mapp  2012; 33: 2072– 80. Google Scholar CrossRef Search ADS PubMed  Weber KA, Sentis AI, Bernadel-Huey ON, Chen Y, Wang X, Parrish TB, et al.   Thermal stimulation alters cervical spinal cord functional connectivity in humans. Neuroscience  2018; 369: 40– 50. Google Scholar CrossRef Search ADS PubMed  Wheeler-Kingshott CA, Stroman PW, Schwab JM, Bacon M, Bosma R, Brooks J, et al.   The current state-of-the-art of spinal cord imaging: applications. Neuroimage  2014; 84: 1082– 93. Google Scholar CrossRef Search ADS PubMed  Winder AT, Echagarruga C, Zhang Q, Drew PJ. Weak correlations between hemodynamic signals and ongoing neural activity during the resting state. Nat Neurosci  2017; 20: 1761– 9. Google Scholar CrossRef Search ADS PubMed  Wu TL, Wang F, Mishra A, Wilson GHIII, Byun N, Chen LM, et al.   Resting-state functional connectivity in the rat cervical spinal cord at 9.4 T. Magn Reson Med  2018; 79: 2773– 83. Google Scholar CrossRef Search ADS PubMed  Xu Z, Conrad BN, Baucom RB, Smith SA, Poulose BK, Landman BA. Abdomen and spinal cord segmentation with augmented active shape models. J Med Imaging  2016; 3: 36002. Google Scholar CrossRef Search ADS   Zhou F, Zhuang Y, Gong H, Wang B, Wang X, Chen Q, et al.   Altered inter-subregion connectivity of the default mode network in relapsing remitting multiple sclerosis: a functional and structural connectivity study. PLoS One  2014; 9: e101198. Google Scholar CrossRef Search ADS PubMed  © The Author(s) (2018). Published by Oxford University Press on behalf of the Guarantors of Brain. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Brain Oxford University Press

Multiple sclerosis lesions affect intrinsic functional connectivity of the spinal cord

Loading next page...
 
/lp/ou_press/multiple-sclerosis-lesions-affect-intrinsic-functional-connectivity-of-fhfJ6Gm9bD
Publisher
Oxford University Press
Copyright
© The Author(s) (2018). Published by Oxford University Press on behalf of the Guarantors of Brain. All rights reserved. For permissions, please email: journals.permissions@oup.com
ISSN
0006-8950
eISSN
1460-2156
D.O.I.
10.1093/brain/awy083
Publisher site
See Article on Publisher Site

Abstract

Abstract Patients with multiple sclerosis present with focal lesions throughout the spinal cord. There is a clinical need for non-invasive measurements of spinal cord activity and functional organization in multiple sclerosis, given the cord’s critical role in the disease. Recent reports of spontaneous blood oxygenation level-dependent fluctuations in the spinal cord using functional MRI suggest that, like the brain, cord activity at rest is organized into distinct, synchronized functional networks among grey matter regions, likely related to motor and sensory systems. Previous studies looking at stimulus-evoked activity in the spinal cord of patients with multiple sclerosis have demonstrated increased levels of activation as well as a more bilateral distribution of activity compared to controls. Functional connectivity studies of brain networks in multiple sclerosis have revealed widespread alterations, which may take on a dynamic trajectory over the course of the disease, with compensatory increases in connectivity followed by decreases associated with structural damage. We build upon this literature by examining functional connectivity in the spinal cord of patients with multiple sclerosis. Using ultra-high field 7 T imaging along with processing strategies for robust spinal cord functional MRI and lesion identification, the present study assessed functional connectivity within cervical cord grey matter of patients with relapsing-remitting multiple sclerosis (n = 22) compared to a large sample of healthy controls (n = 56). Patient anatomical images were rated for lesions by three independent raters, with consensus ratings revealing 19 of 22 patients presented with lesions somewhere in the imaged volume. Linear mixed models were used to assess effects of lesion location on functional connectivity. Analysis in control subjects demonstrated a robust pattern of connectivity among ventral grey matter regions as well as a distinct network among dorsal regions. A gender effect was also observed in controls whereby females demonstrated higher ventral network connectivity. Wilcoxon rank-sum tests detected no differences in average connectivity or power of low frequency fluctuations in patients compared to controls. The presence of lesions was, however, associated with local alterations in connectivity with differential effects depending on columnar location. The patient results suggest that spinal cord functional networks are generally intact in relapsing-remitting multiple sclerosis but that lesions are associated with focal abnormalities in intrinsic connectivity. These findings are discussed in light of the current literature on spinal cord functional MRI and the potential neurological underpinnings. multiple sclerosis, spinal cord, lesions, functional connectivity, resting state functional magnetic resonance imaging Introduction Clinical background The spinal cord is a critical structure in health and disease. Damage to the cord can lead to dysfunction of motor and/or sensory systems, as occurs in cases of trauma involving spinal cord injury, as well as in various musculoskeletal and neurodegenerative diseases such as degenerative disc disease and multiple sclerosis. Spinal cord dysfunction may also be implicated in the absence of gross structural damage, such as in chronic pain syndrome. Given the central role of spinal cord dysfunction in these patient populations, there is an important (but largely unmet) clinical need for non-invasive measurements of cord activity and functional organization (Wheeler-Kingshott et al., 2014). Functional connectivity in the spinal cord We recently reported the presence of resting state blood oxygenation level-dependent (BOLD) correlations in the human spinal cord using ultra-high field MRI, as well as evidence that these effects are reproducible within and across human subjects (Barry et al., 2014, 2016). Such observations have since been corroborated by several research groups using alternative acquisition, processing, and statistical methods (Kong et al., 2014; Liu et al., 2016; San Emeterio Nateras et al., 2016; Eippert et al., 2017b; Harita and Stroman, 2017; Weber et al., 2018), as well as in non-human primate and rat spinal cord (Chen et al., 2015; Wu et al., 2018). These studies have demonstrated that spontaneous BOLD signal fluctuations in the spinal cord are primarily characterized by ventral motor and dorsal sensory networks. The organization and spatial distribution of this activity at rest is also consistent with task-based functional MRI studies (Maieron et al., 2007; Cadotte et al., 2012; Stroman et al., 2012). Eippert et al. (2017b) looked at resting state correlations in cervical cord grey matter at 3 T using multiple variations on their preprocessing and analysis pipelines, finding that the ventral-ventral and dorsal-dorsal horn networks were robust to variations in processing procedures, successfully replicating our findings at 7 T. However, it is important to examine how spontaneous spinal cord activity and functional connectivity is altered in disease states, particularly diseases that result in cord dysfunction. Overall goals In the current study, we (i) examined resting state BOLD fluctuations in the cervical spinal cord grey matter in patients with relapsing-remitting multiple sclerosis compared to a large sample of healthy controls; (ii) assessed the relationship of spinal cord connectivity measures with standard clinical measures of disability; and (iii) assessed the effects of lesions on region-to-region connectivity both locally and distant in relation to ascending/descending white matter tracts. The present study provides the first assessment of functional MRI-derived functional connectivity in the cord of a diseased population. Structural and functional organization of the spinal cord The spinal cord is responsible for the transmission and conduction of electrical potentials to and from the brain, the primary interface of the central and peripheral nervous systems. It is well organized and is composed of grey and white matter regions. The spinal cord grey matter is centrally located surrounded by white matter and is characterized by its ‘butterfly’ shape. Spinal cord white matter contains dense bundles of myelinated axons running in the rostrocaudal direction, separated into the dorsal, lateral, and ventral columns. Spinal cord grey matter can be divided into functional subunits corresponding to motor and sensory inputs, including anterior (ventral) and posterior (dorsal) horns (Cohen-Adad and Wheeler-Kingshott, 2014). Motor network Descending spinal pathways originating primarily in motor cortex are responsible for voluntary motor control. Projections from contralateral cortical sites decussate (cross hemispheres) in the medulla of the brainstem and travel down the cord along ventral and lateral corticospinal tracts. Corticospinal tract axons synapse onto neurons in the ventral grey matter of the spinal cord, from which efferent fibres exit the ventral roots of the cord and innervate specific muscle groups (Purves et al., 2012; Kandel et al., 2013). Sensory network Ascending spinal pathways carry somatosensory information up to the brain from a myriad of peripheral afferents, including tactile, temperature, and pain receptors in the skin as well as from muscles, joints, and visceral organs. Mechanoreceptors and proprioceptors project onto the dorsal column medial lemniscus (DCML) pathway and are responsible for encoding information about fine touch, vibration, pressure, and position. These receptors send their first-order projections from ipsilateral dorsal root ganglia (location of sensory receptor cell bodies) to the spinal cord where the axons travel along ipsilateral dorsal column tracts and ultimately synapse onto grey matter neurons in the medulla. Alternatively, nociceptors and thermoreceptors, which encode information about pain, temperature, and coarse touch, project onto the spinothalamic pathway, which involves anterolateral white matter tracts. For spinothalamic afferents, first-order axons from ipsilateral dorsal root ganglia enter the spinal cord and then synapse onto grey matter neurons in the ipsilateral dorsal horn of the cord. Interneurons and central pattern generators Local connections exist within spinal cord segments that may demonstrate distinct, synchronized activity. The existence of central pattern generators (CPGs) in the spinal cord was observed over 100 years ago (Brown, 1911; McCrea and Rybak, 2008). In the spinal cord, CPGs refer to neural computations of locomotor actions and reflexes, which are performed within the spinal cord and without specific input from higher-order brain areas. These computations are mediated by interneurons that connect local grey matter regions, allowing, for instance, sensory input to locally influence ventral motor neuron output. Challenges in spinal cord functional MRI While functional MRI studies of the brain have increased dramatically over the past 25 years, the detection of robust BOLD signal in the spinal cord has proved challenging (Giove et al., 2004; Stroman et al., 2014). The cervical spinal cord is ∼1–1.5 cm in diameter and grey matter butterfly horns only a couple of millimetres across, and the typical spatial resolution of 2–3 mm used in brain functional MRI is inadequate for spinal cord functional MRI. Higher in-plane resolution is needed to accurately distinguish spatial properties of BOLD signals across the spinal cord and to reduce partial volume effects from nearby tissue such as white matter and CSF. Additionally, susceptibility distortions arising from single-shot echo-planar imaging (EPI), used in brain functional MRI sequences are amplified in spinal cord due to the close proximity to surrounding bones and intervertebral discs and the oesophagus and lungs. A further complication in measuring the spinal cord BOLD signal is an increased influence of physiological noise compared to the brain (Eippert et al., 2017a). These noise sources include: cord motion, respiration and swallowing, transient B0 inhomogeneity, flow and CSF pulsation. We have developed improved acquisition and processing strategies for robust spinal cord functional MRI (Barry et al., 2011, 2013, 2014, 2016, 2018; Conrad et al., 2017), which are used in this study and outlined in the ‘Materials and methods’ section. Functional connectivity of sensorimotor networks in multiple sclerosis Given the spinal cord’s high degree of axonal connections to sensorimotor cortices, as well as reports of synchronous task-based BOLD activity in spinal cord with these brain regions, it is likely that spinal cord resting state BOLD activity is associated with sensorimotor network activity. Studies looking at functional connectivity in multiple sclerosis in the brain have demonstrated a number of abnormalities compared to healthy controls, with reports of both increased and decreased correlation among networks. Further evidence suggests that these abnormalities may be dynamic over time, with compensatory increases in connectivity occurring early in the disease and decreased connectivity with the accumulation of structural damage (Roosendaal et al., 2010; Faivre et al., 2012). It is possible that analogous changes occur in the sensorimotor network with improvement of motor deficits in response to treatment, for example. Such improvements may manifest in and/or involve spinal cord circuits. While more work is needed, the growing body of evidence supporting widespread network alterations in multiple sclerosis suggests resting state functional MRI measures may provide useful indicators of disease pathology and progression (Rocca et al., 2012b, 2015; Filippi and Rocca, 2013; Sacco et al., 2013). Complementary to studies of the brain, the investigation of intrinsic BOLD fluctuations and functional connectivity in the spinal cord may thus provide a more comprehensive understanding of nervous system dysfunction in multiple sclerosis. Furthermore, measures of functional connectivity in the cord may provide surrogate markers of sensorimotor deficits and aid in disease monitoring or prognosis. Functional imaging of the spinal cord in multiple sclerosis While we present new findings in spinal cord resting state functional MRI in patients with multiple sclerosis, several previous studies have looked at task-based functional activation of the cervical cord in multiple sclerosis in response to tactile and proprioceptive stimulation of the hand or wrist (Agosta et al., 2008a, b, 2009; Valsasina et al., 2010, 2012). This body of work has shown that patients with multiple sclerosis demonstrate increased activation of the cervical cord compared to control subjects, as well as a more bilateral response in the context of ipsilateral stimulation. We may expect over-recruitment of bilateral networks in a task-paradigm to manifest as increased connectivity in a resting state, and therefore predicted that patients with multiple sclerosis would demonstrate an overall increase in connectivity between bilateral region pairs in comparison to controls. Importantly, abnormalities detected via spinal cord functional MRI may relate to clinical measures in multiple sclerosis. One study looking at the relationship of fatigue in patients with multiple sclerosis with cord activity in response to palm stimulation found increased activation in patients with low fatigue scores as compared to controls (Rocca et al., 2012a). Interestingly, a second group of patients with high fatigue scores showed reduced activation compared to those with low fatigue, and instead demonstrated a more diffuse pattern of activity compared to controls. Similar to findings in the brain in multiple sclerosis, it is possible that abnormal functional activity in the cord may follow a dynamic trajectory of early compensation followed by decreased recruitment in later stages concomitant with increasing disability and structural damage. In addition to assessing global abnormalities in spinal cord functional connectivity in multiple sclerosis, a primary aim of the present study was to assess the effects of spinal cord lesion location in multiple sclerosis on functional connectivity among grey matter regions of the cervical cord. Dula et al. (2016) recently showed that 50% more lesions were detectable with 7 T imaging compared to 3 T using similar scan times in the same group of patients. The present study takes advantage of this increased sensitivity to lesions through the use of 7 T imaging. In regards to prior work on spinal cord connectivity and sites of local damage such as with multiple sclerosis lesions, a study by Chen et al. (2015) looked at functional connectivity in the cervical cord of anaesthetized squirrel monkeys before and after introducing a unilateral dorsal lesion in the context of a spinal cord injury model. Connectivity was reduced within the slice of the lesion as well as slices below the lesion involving ipsilateral ventral and dorsal horns, whereas connectivity above the lesion and in contralateral horns was largely unaffected. These results demonstrate that focal lesions in the cord may have both local (intrasegmental) as well as distal (intersegmental) effects on spontaneous fluctuations in spinal grey matter activity. Recovery of disrupted connectivity towards pre-injury levels was observed in a time period of 7–24 weeks after lesioning. This suggests that reorganization of spinal circuitry occurs in response to damage and that monitoring of resting state connectivity may serve as a useful clinical indicator of functional integrity. It is important to note that lesions in the spinal cord are less well studied than those in the brain, largely due to limitations in the spatial resolution of MRI at clinical field strengths. Some work has shown that specific lesion locations, such as within the corticospinal tracts and motor network in general, are more predictive of disease progression (Bodini et al., 2011), providing further motivation for assessments of spinal cord lesions and their impact on functional integrity. Materials and methods Subjects and image acquisition All subjects provided signed, informed consent prior to inclusion in the study in accordance with the local institutional review board. Patients with relapsing-remitting multiple sclerosis (n = 22) were recruited in collaboration with the Vanderbilt Multiple Sclerosis Center as part of a multi-modal 7 T MRI study. Healthy volunteers (n = 56) were also recruited and the first resting state run was used for analysis if two were acquired. Subject demographics for our patient and control groups are provided in Table 1. An expanded set of clinical variables and symptoms for each patient is provided in Supplementary Table 1. All imaging was performed on a Philips Achieva 7 T scanner using a cervical spine coil (Nova Medical Inc.) with quadrature transmit and a 16-channel receive array. Foam padding was used in positioning subjects on the coil such that the cervical spine was straight but comfortable. Subjects were instructed to lie still throughout the scan. Physiological monitoring was performed for all subjects and involved recording from a pulse oximeter placed on a finger and respiration bellow placed on the abdomen. Table 1 Subject demographics Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  EDSS = Expanded Disability Status Scale; TSNR = temporal signal-to-noise ratio. Table 1 Subject demographics Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  Group  Age (years)  Gender  Handedness  TSNR  EDSS  Multiple sclerosis duration (years)  25 ft Walk  Controls, n = 56  29.2 ± 9.8  29F/27M  50R/6L  34.8 ± 5.9  –  –  –  Patients, n = 22  41.6 ± 9.1  12F/10M  20R/2L  33.1 ± 8.9  2.3 ± 2.1  9.9 ± 7.6  6.2 ± 2.7  EDSS = Expanded Disability Status Scale; TSNR = temporal signal-to-noise ratio. For each subject, a high resolution sagittal T2-weighted image was initially acquired and used to plan axial imaging such that the axial slice stack was centred on the C3-C4 vertebral disc and perpendicular to the spinal cord. Curvature of the cervical spine was sometimes not fully mitigated by padding and thus an orthogonal slice placement was not possible. In such cases, slice angulation was planned such that the middle slices were perpendicular to the cord at the C3-C4 vertebral level. Prior to functional imaging, a T2*-weighted anatomical volume was acquired with the same axial slice geometry as the functional scan and the following parameters: field of view = 160 × 160 mm2, 12 4-mm slices (centred on the C3/C4 junction), nominal voxel size = 0.6 × 0.6 × 4 mm3, interpolated voxel size = 0.31 × 0.31 × 4 mm3, repetition time = 303 ms, echo time = 8.2 ms, flip angle = 25°, sensitivity encoding (SENSE) (Pruessmann et al., 1999) reduction factor = 2.0 (right-left), number of acquisitions = 8, total acquisition time = 5 min and 22 s. The anatomical images provided high resolution detail of the cord, including grey matter, white matter, and CSF boundaries, and were used to delineate structures of interest for functional data processing. Additionally, the T2*-weighted anatomical images were used for lesion ratings in patients. For functional imaging, a T2*-weighted resting state functional MRI scan with 150 dynamics was acquired, using a 3D multishot gradient echo sequence with the following parameters: field of view = 160 × 160 mm2, 12 4-mm slices, voxel size = 0.91 × 0.91 × 4 mm3, repetition time = 17 ms, echo time = 8.0 ms, flip angle = 15°, echo train length (number of k-space lines acquired per shot) = 9, SENSE reduction factor = 1.56 (anterior-posterior), volume acquisition time = 3.34 s (278 ms/slice), number of volumes = 150 (after 10 ‘dummy’ scans), total scan time ∼9 min. Preprocessing was the same for each subject and followed a pipeline previously described in detail (Barry et al., 2016), involving motion correction, physiological signal de-noising, white matter and CSF signal regression, and bandpass filtering from 0.01–0.13 Hz (Cox, 1996; Glover et al., 2000). See the online Supplementary material for more information on these procedures. Region of interest correlation analysis The grey matter masks were divided into four quadrants (regions of interest) (Fig. 1) representing the ventral and dorsal horns, excluding the central grey matter commissure and central canal. The grey matter regions of interest were then morphologically eroded to remove the outermost layer of interpolated voxels, minimizing partial volume effects from white matter. The mean signal time course across all voxels within each region of interest was extracted. Pearson correlation coefficients were computed for each of the six region of interest pairs within each slice. The correlation coefficients were transformed to z-scores using a Fisher transformation: z = tanh−1(r)(dof − 3)1/2, where dof is the estimated degrees of freedom after correction for first-order autocorrelation (Rogers and Gore, 2008) and z represents the distance from zero correlation in units of standard deviation for each region of interest pair. Note that while we have explored the use of partial correlations in previous work (i.e. including as nuisance variables the time series from the regions of interest opposite to the pair being correlated), we chose to analyse these data with full correlations to simplify interpretation of the results. Figure 1 View largeDownload slide Imaging and resting state connectivity.Top left: Example slice placement as centred over C3-C4 vertebral disk, used as landmark for placement in all subjects. Top right: Example high resolution, T2*-weighted anatomical image and close-up; interp. = interpolated voxel size (displayed in figure). Middle right: Example functional MRI-weighted acquisition and schematic of functional imaging showing 150 volumes acquired over time (t), which was every 3.34 s. Functional images acquired using a 3D multishot gradient echo sequence, images shown after registration and interpolation to anatomical space. Middle left: Example seed-to-voxel correlation map (thresholded at Pearson’s r > 0.465, P < 10−7 uncorrected). An example seed voxel in left ventral grey matter (GM) shows correlation with contralateral ventral, ipsilateral dorsal, and contralateral dorsal grey matter. Bottom left: Spinal cord regions of interest defined during manual segmentation and preprocessing. CC = central commissure of grey matter (voxels excluded from analyses); LD = left dorsal grey matter horn; LV = left ventral grey matter horn; RD = right dorsal grey matter horn; RV = right ventral grey matter horn; NS = ‘not-spine’ mask (included all voxels outside of the spinal canal); WM = white matter. Bottom right: Representative functional MRI BOLD signal timecourses in a single slice from LV and RV, and LD and RD, respectively, after bandpass filter. Pearson’s r and associated Z-score for correlation between both region of interest pairs provided. Figure 1 View largeDownload slide Imaging and resting state connectivity.Top left: Example slice placement as centred over C3-C4 vertebral disk, used as landmark for placement in all subjects. Top right: Example high resolution, T2*-weighted anatomical image and close-up; interp. = interpolated voxel size (displayed in figure). Middle right: Example functional MRI-weighted acquisition and schematic of functional imaging showing 150 volumes acquired over time (t), which was every 3.34 s. Functional images acquired using a 3D multishot gradient echo sequence, images shown after registration and interpolation to anatomical space. Middle left: Example seed-to-voxel correlation map (thresholded at Pearson’s r > 0.465, P < 10−7 uncorrected). An example seed voxel in left ventral grey matter (GM) shows correlation with contralateral ventral, ipsilateral dorsal, and contralateral dorsal grey matter. Bottom left: Spinal cord regions of interest defined during manual segmentation and preprocessing. CC = central commissure of grey matter (voxels excluded from analyses); LD = left dorsal grey matter horn; LV = left ventral grey matter horn; RD = right dorsal grey matter horn; RV = right ventral grey matter horn; NS = ‘not-spine’ mask (included all voxels outside of the spinal canal); WM = white matter. Bottom right: Representative functional MRI BOLD signal timecourses in a single slice from LV and RV, and LD and RD, respectively, after bandpass filter. Pearson’s r and associated Z-score for correlation between both region of interest pairs provided. Lesion ratings Lesion ratings were independently conducted by three raters (B.N.C., S.M., and A.B.) who were blind to patient clinical information during the rating process. Control images of the spinal cord were also provided to the raters to determine the normal signal intensity of the axial T2*-weighted images, e.g. hyperintensities present in healthy white matter, in an effort to eliminate false-positive lesion identification in the patient data. For every slice (264 total slices across all patients), raters determined whether a lesion was present in each of the four white matter columns: ventral, dorsal, left lateral, and/or right lateral, yielding a total of 1056 possible lesion locations. Raters also had the option to not rate an entire slice if an artefact interfered with lesion demarcation. After independent ratings were completed, a consensus meeting was held in which slices were reassessed in the presence of all three raters. Discrepancies in the original ratings were reviewed and a final rating was agreed upon in each instance. Statistical analyses Group differences of connectivity values in patients compared to controls were assessed with two-tailed Wilcoxon rank-sum tests, performed using the mean connectivity value across all 12 slices for each subject. To assess the effect of lesions on slicewise connectivity values, we used linear mixed effects models in which the region of interest pair connectivity z-scores were the dependent variables. A random effect of subject was included (since connectivity values were generally higher in some subjects compared to others), and then fixed effects of subject group (control versus patient), age, gender, median slicewise temporal signal-to-noise ratio (TSNR) in grey matter, and the four binary lesion location variables (controls with all zeros). Additional analyses (not shown) including handedness as a covariate showed no significant main effects of handedness and did not appreciably affect the results presented herein. Results Connectivity in controls Connectivity results for the six region of interest pairings (where V = ventral horn, D = dorsal horn, R = right, L = left) across all slices in healthy controls (672 total slices, grey dots) are shown in Fig. 2A. Ventral-ventral (V-V) connectivity was the highest of all pairings, followed by dorsal-dorsal (D-D), and then the ipsilateral pairings. The diagonal ventral-dorsal pairings demonstrated the lowest connectivity overall. Connectivity across vertebral level (V-V and D-D) is plotted in Fig. 2B and C. There were no significant differences in V-V connectivity across vertebral levels. D-D connectivity for slices within C5 were significantly higher than the other levels obtained in our slice stack using two-tailed Wilcoxon rank-sum tests. In Fig. 2D, Z-scores for the V-V pairing in controls are grouped by gender and plotted for each subject, ordered according to subject median values. The dotted lines indicate the median z-score across all slices for each gender. Females had significantly higher V-V functional connectivity compared to males (for more on this gender effect, see the Supplementary material). Additionally, Fig. 2D demonstrates the subject-level dependence of connectivity values, motivating the linear mixed effects modelling undertaken in subsequent analyses where a random effect of subject was taken into account. Figure 2 View largeDownload slide Functional connectivity in controls and multiple sclerosis patients. (A) Box/dot plots of z-scores (i.e. connectivity strength) from every slice in controls (56 subjects × 12 slices = 672 total) for each region of interest pairing. Z-score = tanh−1(r) (dof − 3)1/2, where r is the Pearson correlation coefficient between mean timeseries from each region of interest and dof is the estimated degrees of freedom after correction for first-order autocorrelation. The mean z-score across all slices is indicated by a white midline, with black indicating the standard error of the mean and boxes indicating the standard deviation. (B and C) Box/dot plots showing z-scores for V-V and D-D across vertebral levels. Boxplots indicators are the same as in A. Two-tailed Wilcoxon rank-sum tests were performed, *P < 0.05, **P < 0.01. (D) Boxplots of z-scores for V-V in each of the 56 control subjects (12 slices/z-scores per subject), grouped by gender and ordered according to median z-score (lowest to highest within the group). The median z-score for each subject is indicated by a dark midline and boxes indicate the 25–75% quartile ranges. Whiskers extend to the most extreme z-scores not considered outliers and outliers are indicated with a black dot. Red and orange dotted lines indicate the median z-score across all males and females, respectively. (E) Group average power, after bandpass filtering. Each subject’s mean power estimated across all grey matter voxels. Wilcoxon rank-sum tests detected no significant differences between groups after division of the spectrum into 10 frequency bins. (F) Boxplots of mean z-scores for each region of interest pairing in controls (n = 56) and patients (n = 22). Subjects contributed a single mean z-score across 12 slices for each region of interest pair. Boxplots indicators are the same as previously described for D. Two-tailed Wilcoxon rank-sum tests detected no significant differences between control and patients; however, there was a trend towards decreased D-D as well as RV-LD connectivity in patients (P = 0.097 and P = 0.095, respectively). Figure 2 View largeDownload slide Functional connectivity in controls and multiple sclerosis patients. (A) Box/dot plots of z-scores (i.e. connectivity strength) from every slice in controls (56 subjects × 12 slices = 672 total) for each region of interest pairing. Z-score = tanh−1(r) (dof − 3)1/2, where r is the Pearson correlation coefficient between mean timeseries from each region of interest and dof is the estimated degrees of freedom after correction for first-order autocorrelation. The mean z-score across all slices is indicated by a white midline, with black indicating the standard error of the mean and boxes indicating the standard deviation. (B and C) Box/dot plots showing z-scores for V-V and D-D across vertebral levels. Boxplots indicators are the same as in A. Two-tailed Wilcoxon rank-sum tests were performed, *P < 0.05, **P < 0.01. (D) Boxplots of z-scores for V-V in each of the 56 control subjects (12 slices/z-scores per subject), grouped by gender and ordered according to median z-score (lowest to highest within the group). The median z-score for each subject is indicated by a dark midline and boxes indicate the 25–75% quartile ranges. Whiskers extend to the most extreme z-scores not considered outliers and outliers are indicated with a black dot. Red and orange dotted lines indicate the median z-score across all males and females, respectively. (E) Group average power, after bandpass filtering. Each subject’s mean power estimated across all grey matter voxels. Wilcoxon rank-sum tests detected no significant differences between groups after division of the spectrum into 10 frequency bins. (F) Boxplots of mean z-scores for each region of interest pairing in controls (n = 56) and patients (n = 22). Subjects contributed a single mean z-score across 12 slices for each region of interest pair. Boxplots indicators are the same as previously described for D. Two-tailed Wilcoxon rank-sum tests detected no significant differences between control and patients; however, there was a trend towards decreased D-D as well as RV-LD connectivity in patients (P = 0.097 and P = 0.095, respectively). Power spectra in patients and controls In Fig. 2E, mean power spectra for all grey matter voxels were estimated for each subject and averaged across the group for both controls (red/orange) and patients (blue), as well as an age- and gender-matched subset of the control group (12 females/10 males, mean age = 38.2 ± 10.2 years) (cyan). For this analysis, a 251-point fast Fourier transform was performed for each region of interest timecourse, yielding 100 points within the frequency range of interest (i.e. range used in bandpass filtering, 0.01–0.13 Hz). To test for differences within specific frequency bands, the range of interest was divided into 10 frequency bins with 10 points per subject. The values within each bin were averaged for each subject and Wilcoxon rank-sum tests were performed to look for group differences. No frequency bin was significantly different between the patients and either of the control groups. The qualitative increase in spectral power observed in the patient group in Fig. 2E is thus likely reflective of outliers rather than a true group difference. The spike in power observed in all groups at approximately 0.075 Hz is characteristic of a low frequency oscillation observed in brain resting state studies (Kiviniemi et al., 2005; Tong et al., 2011). It has been suggested that these spontaneous fluctuations in BOLD signal are reflective of vasomotor oscillations in arterial blood pressure/flow, e.g. Mayer-waves or C-waves (Katura et al., 2006). Fluctuations around the observed 0.075 Hz peak survived in grey matter after white matter denoising (Supplementary Fig. 1) suggesting this frequency range is related to grey matter neuronal activity. However, our white matter denoising approach does not constitute a complete suppression of signals from white matter (and surrounding spinal cord vasculature), and thus we cannot rule out the possibility that the peak at 0.075 Hz represents residual physiological noise, such as from vasomotor oscillations. There is ongoing debate about the neuronal origin of resting state BOLD fluctuations in the brain, with a recent paper suggesting that spontaneous vasomotion in the range of ∼0.1 Hz tracks the envelope of gamma-band electrical activity in the cortex (Mateo et al., 2017), while another suggests that resting state BOLD fluctuations are influenced by unrelated haemodynamic processes and only weakly correlate with neural activity (Winder et al., 2017). Future research is required to further investigate the source of spinal cord BOLD signals and their potential neuronal origins, such as via intraspinal electrical recording. Further analyses assessing the contribution of narrower frequency bands to our functional connectivity metrics are provided in Supplementary Fig. 2. Figure 3 View largeDownload slide Lesion ratings. Patient lesion ratings were conducted by three independent raters and involved determination for each slice whether a lesion was present in ventral, dorsal, left lateral, and/or right lateral column of the spinal cord white matter. Final consensus ratings were determined and are depicted for each subject. Example slices, which were rated as containing a lesion, are provided. The straightened spinal cord and representative cervical slice images to the right were adapted from recently developed spinal cord templates included in the Spinal Cord Toolbox software package (Fonov et al., 2014; De Leener et al., 2017, 2018). Figure 3 View largeDownload slide Lesion ratings. Patient lesion ratings were conducted by three independent raters and involved determination for each slice whether a lesion was present in ventral, dorsal, left lateral, and/or right lateral column of the spinal cord white matter. Final consensus ratings were determined and are depicted for each subject. Example slices, which were rated as containing a lesion, are provided. The straightened spinal cord and representative cervical slice images to the right were adapted from recently developed spinal cord templates included in the Spinal Cord Toolbox software package (Fonov et al., 2014; De Leener et al., 2017, 2018). Comparison of functional connectivity in patients to controls To look for differences in functional connectivity between patient and control groups, the mean z-score across 12 slices for each subject/region of interest pairing was first computed (Fig. 2F). Tests of difference were conducted on the mean values using two-tailed Wilcoxon rank-sum tests for each region of interest pairing. There were no significant differences in z-scores between subject groups for any of the region of interest pairs, although a trend towards decreased D-D as well as RV-LD connectivity was observed in patients (P = 0.097 and P = 0.095, respectively). While the rank-sum tests served as a first pass for looking at group differences in the raw connectivity values, they did not control for other variables such as the significant difference in age between the patient and control groups or potential gender effects, as is demonstrated in Fig. 2D. To further assess group differences while controlling for other variables, linear models were fit using multiple regression. Four predictor variables were included in each model, group (Control/Patient), age, gender (M/F), and TSNR (calculated after white matter denoising but prior to bandpass filtering). The main effect of each variable in predicting the average region of interest pair z-score over 12 slices is depicted in Supplementary Fig. 4 along with the confidence intervals for these values. None of the region of interest pairings demonstrated a significant main effect of group. One challenge in interpreting these results is that our patient sample was significantly older, potentially conflating effects of disease and normal ageing, for instance in the D-D connectivity. Finally, we assessed correlations between clinical measures of disability and average connectivity for each region of interest pairing in the patient group. We found that no disability measure [i.e. Expanded Disability Status Scale (EDSS), disease duration, and timed 25-foot walk test] was significantly directly correlated with mean connectivity strength in patients (Supplementary Fig. 3). When controlling for nuisance variables (i.e. age, gender, and TSNR), however, we found modest positive associations of EDSS and 25-foot walk times with D-D connectivity, as well a positive association between 25-foot walk times and V-V connectivity (Supplementary Fig. 5). These findings indicate that increasing disability is related to increased bilateral connectivity, and it is possible that this reflects a compensatory mechanism among bilateral spinal cord networks. However, these observations should be considered preliminary and more work with a larger sample size is necessary to confirm this finding. Lesion ratings analysis To assess the level of agreement between the three independent raters, we looked at the per cent agreement across all possible lesions locations (i.e. lesion present or not) in which both raters provided a rating (note that raters had the option to not rate an entire slice due to artefact). Between each pair of raters, the per cent agreement was 90.7% (1028 total locations), 88.2% (952 total locations), and 91.4% (940 total locations), respectively. Following this analysis, a consensus meeting was held in which the discrepancies were reviewed and final ratings were agreed upon. In five slices (out of 264, ∼2%), a rating was not assigned because of artefacts. Of the possible 1036 locations considered for the final analysis, 265 locations were rated as containing a lesion. Overall, 19 of 22 patients presented with a lesion somewhere within the imaged volume. The final consensus lesion ratings were taken as the ground truth for subsequent analyses and are visually depicted in Fig. 3. Local effects of lesions on functional connectivity To assess the effect of lesions on connectivity within a slice (local), linear mixed effects models were fit individually for each of the region of interest pairings, with the slicewise connectivity z-score as the outcome variable, a random effect of subject, and the remaining variables as fixed effects. Control data were included in these models, with lesion values at zero. The relevant statistical results for the six models are provided in Table 2, and the full statistical data including parameter estimates for nuisance variables are provided in Supplementary Table 2. Effect sizes and polarity of the estimate (whether positive or negative) for each lesion location are visually depicted in Fig. 4. Ventral column lesions demonstrated a significant negative association with local V-V connectivity such that the presence of a ventral lesion was associated with decreased V-V connectivity after controlling for all other effects, including the effects of lesions in other locations. In the D-D model, the presence of a left lateral lesion demonstrated a significant positive association with D-D connectivity. Additionally, higher RV-RD connectivity was associated with the presence of right lateral lesions and higher LV-RD connectivity was associated with dorsal column lesions. Table 2 Local effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Table 2 Local effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  Name  Estimate  P-value  Lower  Upper  V-V          L lesion  0.08  0.857  −0.79  0.95  R lesion  −0.30  0.531  −1.25  0.65  D lesion  −0.47  0.329  −1.41  0.47  V lesion  −1.17  0.043*  −2.30  −0.04  D-D          L lesion  0.88  0.043*  0.03  1.74  R lesion  0.66  0.171  −0.28  1.60  D lesion  0.15  0.746  −0.76  1.07  V lesion  −0.37  0.516  −1.49  0.75  LV-LD          L lesion  0.38  0.289  −0.32  1.08  R lesion  0.08  0.841  −0.70  0.85  D lesion  −0.19  0.604  −0.91  0.53  V lesion  −0.48  0.310  −1.40  0.45  RV-RD          L lesion  −0.40  0.262  −1.10  0.30  R lesion  0.87  0.028*  0.09  1.64  D lesion  0.71  0.053  −0.01  1.43  V lesion  −0.02  0.961  −0.95  0.90  LV-RD          L lesion  0.46  0.200  −0.25  1.18  R lesion  0.07  0.863  −0.72  0.86  D lesion  0.78  0.036*  0.05  1.52  V lesion  −0.03  0.954  −0.97  0.91  RV-LD          L lesion  0.19  0.595  −0.51  0.89  R lesion  −0.02  0.964  −0.80  0.76  D lesion  0.54  0.144  −0.18  1.26  V lesion  −0.55  0.242  −1.48  0.37  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Figure 4 View largeDownload slide Local effects of lesions on functional connectivity. Visual representation of effect sizes for lesion presence (ventral, dorsal, left lateral, and right lateral column) in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 2. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue circles indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the column, which they represent. All patient (n = 22) and control (n = 56) data were included in the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for lesion values. Full statistical results are provided in Supplementary Table 2. Note that all lesion variables were included in the same model, such that the estimates provided are the independent effect of a lesion within the column on connectivity after accounting for all other effects, including lesions in other columns. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes in table indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. FC = functional connectivity. Figure 4 View largeDownload slide Local effects of lesions on functional connectivity. Visual representation of effect sizes for lesion presence (ventral, dorsal, left lateral, and right lateral column) in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 2. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue circles indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the column, which they represent. All patient (n = 22) and control (n = 56) data were included in the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for lesion values. Full statistical results are provided in Supplementary Table 2. Note that all lesion variables were included in the same model, such that the estimates provided are the independent effect of a lesion within the column on connectivity after accounting for all other effects, including lesions in other columns. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes in table indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. FC = functional connectivity. Distal effects of lesions on functional connectivity To assess the effect of lesions on connectivity above and below the lesion, we first classified lesions as occurring in ascending and/or descending white matter tracts (see schematic in Fig. 5). Dorsal column lesions were classified as ascending tract lesions. Ventral column lesions were classified as descending tract lesions. Since both ascending and descending tracts are found in the lateral columns, these lesions were counted for both classes. Our anatomical imaging was not of sufficient resolution to confidently distinguish lesions within specific lateral column tracts. Furthermore, many of the lateral lesions encompassed large areas of the column, thus likely affecting both ascending and descending tracts (see lesion examples in Fig. 3). Slices were then assigned a value of 1 if they occurred either above an ascending tract lesion or below a descending tract lesion. Linear mixed effects models were then individually fit for each of the region of interest pairs and the relevant results are provided in Table 3. The full statistical results including parameter estimates for nuisance variables are provided in Supplementary Table 3. A visual depiction of the effect sizes for lesion location are displayed in Fig. 5. Though not significant at P < 0.05, a modest negative effect was observed for V-V connectivity in slices below a descending tract lesion after controlling for other variables (P = 0.131). LD-RV connectivity in slices below a descending tract lesion was significantly decreased (P = 0.011). Within the same model, an increase in LD-RV connectivity in slices above an ascending tract lesion approached significance (P = 0.057). Table 3 Distal effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Table 3 Distal effects of lesions on functional connectivity Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  Name  Estimate  P-value  Lower  Upper  V-V          Below  −0.69  0.131  −1.59  0.21  Above  −0.32  0.457  −1.16  0.52  D-D          Below  −0.34  0.435  −1.21  0.52  Above  −0.16  0.707  −0.98  0.66  LV-LD          Below  0.19  0.577  −0.47  0.84  Above  −0.11  0.739  −0.76  0.54  RV-RD          Below  0.13  0.707  −0.53  0.78  Above  0.32  0.331  −0.33  0.98  LV-RD          Below  0.25  0.462  −0.42  0.92  Above  0.40  0.245  −0.27  1.06  LD-RV          Below  −0.84  0.011*  −1.50  −0.19  Above  0.64  0.057  −0.02  1.29  *P < 0.05. D = dorsal; L = left; R = right; V = ventral. Figure 5 View largeDownload slide Distal effects of lesions on functional connectivity. Visual representation of effect sizes for distal presence of lesion in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 3. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue arrows indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the regions which they represent. Slices were rated as containing a descending tract lesion if a lesion was present in at least one of the ventral, left lateral, or right lateral tracts; and a descending tract lesion if a lesion was present in at least one of the dorsal, left lateral, or right lateral tracts. Based on this classification, two independent variables were constructed in which slices were assigned a binary value indicating whether they did (1) or did not (0) occur (i) above (superior to) an ascending tract lesion (‘Above’ variable in Table 3, upward facing arrows in figure); and (ii) below (inferior to) a descending tract lesion (‘Below’ variable in Table, downward facing arrows in figure). Slices could receive only binary classifications of 1 or 0 for the Above and Below variables and thus were not additive, for instance in cases where a slice was above multiple ascending tract lesions. All patient (n = 22) and control (n = 56) data were included the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for Above and Below values. Full statistical results are provided in Supplementary Table 3. Note that the Above and Below variables were included in the same model, such that the estimates provided are the independent effect of a slice occurring above an ascending tract lesion after accounting for all other effects, including the slice occurring below a descending tract lesion. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. Figure 5 View largeDownload slide Distal effects of lesions on functional connectivity. Visual representation of effect sizes for distal presence of lesion in relation to region of interest pair z-score. Effect size estimates came from independent linear mixed-effects models of slicewise region of interest pair connectivity strength (z-scores), and are reported in Table 3. In each model, a random effect of subject was included, along with fixed effects of group, age, gender, and median slicewise temporal signal-to-noise ratio in grey matter. Red and blue arrows indicate positive and negative effect estimates, respectively, are proportional to the magnitude of the effect size, and are located in the regions which they represent. Slices were rated as containing a descending tract lesion if a lesion was present in at least one of the ventral, left lateral, or right lateral tracts; and a descending tract lesion if a lesion was present in at least one of the dorsal, left lateral, or right lateral tracts. Based on this classification, two independent variables were constructed in which slices were assigned a binary value indicating whether they did (1) or did not (0) occur (i) above (superior to) an ascending tract lesion (‘Above’ variable in Table 3, upward facing arrows in figure); and (ii) below (inferior to) a descending tract lesion (‘Below’ variable in Table, downward facing arrows in figure). Slices could receive only binary classifications of 1 or 0 for the Above and Below variables and thus were not additive, for instance in cases where a slice was above multiple ascending tract lesions. All patient (n = 22) and control (n = 56) data were included the models, with each subject contributing 12 slicewise values (thus a total of 936 values for each variable, with repeated measures for subject, group, age, and gender). For all slices, controls received zeros for Above and Below values. Full statistical results are provided in Supplementary Table 3. Note that the Above and Below variables were included in the same model, such that the estimates provided are the independent effect of a slice occurring above an ascending tract lesion after accounting for all other effects, including the slice occurring below a descending tract lesion. Green lines/circles indicate the region of interest pair assessed in the model. Red boxes indicate P-values < 0.05. Lower/Upper = lower and upper bounds of confidence interval. Discussion The present study provides the first assessment of functional MRI-derived functional connectivity in the cord of a diseased population. This builds upon recent reports of intrinsic resting state BOLD activity and correlation in the healthy human spinal cord. Spinal cord functional MRI presents unique challenges, such as requiring higher resolution than studies in the brain due to the small structures of interest. We took advantage of ultra-high field 7 T imaging, which provides increased BOLD contrast as well as increased sensitivity in the detection of multiple sclerosis lesions. Our previous work provided a methodological framework for the measurement of functional connectivity in patients, including a 3D-multishot functional MRI sequence, which reduced susceptibility distortions and an optimized preprocessing and denoising pipeline. The sample of healthy controls considered in this analysis represents the largest described in the spinal cord functional MRI literature to date and further supports the existence of distinct and robust functional networks within human cervical cord grey matter. In both patients and controls, the strongest correlations were observed among bilateral ventral grey matter regions, followed by those among dorsal grey matter regions. Importantly, while more research is needed to confirm the neuronal origin of spinal cord BOLD signals, our work suggests spinal cord networks are detectable in the absence of a particular task or stimulation, making the acquisition and translation of these measures more feasible for future clinical application. We found no significant differences in mean connectivity between multiple sclerosis patients and controls at the group level (Fig. 2F), but did detect both decreased and increased connectivity among grey matter regions in relation to local multiple sclerosis lesions (Fig. 4 and Table 2). For instance, we found a significant decrease in V-V connectivity in the presence of ventral column lesions, suggesting that damage to the ventral column white matter may be particularly detrimental to motor network synchronization within the spinal cord. It is possible that ventral column lesions characterize a direct disconnection between intra-segmental ventral grey matter circuits. For other region pairs in which we found increased connectivity near lesion sites (D-D with L lesion, RV-RD with R lesion and trend with D lesion, and LV-RD with D lesion) one hypothesis is that increased connectivity in the presence of lesions and/or axonal loss (Kearney et al., 2015) reflects a compensatory mechanism within local grey matter circuits. One study by Agosta et al. (2008b) found that there was increased bilateral recruitment of grey matter regions in patients with multiple sclerosis in response to stimulation, and that an increased signal change to proprioceptive stimulation was associated with reduced fractional anisotropy in the cord, suggesting that increasing white matter damage is related to a compensatory increase in activity and/or metabolic demand in functioning grey matter. Furthermore, the present findings are in-line with converging evidence of compensatory increases in brain connectivity in relapsing-remitting multiple sclerosis in sensorimotor network as well as other resting state networks, particularly those showing increased connectivity with increased structural damage. For instance, a study using resting functional MRI and diffusion tensor imaging in minimally disabled relapsing-remitting patients compared to controls found that regions of the default mode network were functionally overconnected while structurally disconnected, due to damage in long-range white matter tracts, suggesting a compensatory mechanism of a network, albeit at a more global scale (Zhou et al., 2014). Our findings of increased local connectivity for some region pairs/lesion sites may too reflect a compensatory effect of white matter damage. Though not mutually exclusive, an alternative mechanism has been proposed to explain overactivation and increased bilateral recruitment of cervical cord grey matter in multiple sclerosis, which involves a disruption of inhibitory spinal interneurons, such that there is failure to inhibit a widespread network of primary spinal neurons. Histopathological assessment of multiple sclerosis spinal cord tissue has demonstrated both a reduced number and decreased cell size of interneurons in the upper cervical levels compared to control tissue (Gilmore et al., 2009). An abnormal increase in the intrinsic coactivation (i.e. connectivity) of grey matter in a resting state may also be reflective of interneuron dysfunction. If so, our findings may indicate interneuron effects on intrinsic activity are most pronounced near lesion sites. While further work is necessary to distinguish between these hypotheses, our findings suggest that the functional consequences of multiple sclerosis cord lesions as measured through fluctuations in grey matter BOLD signals, such as from compensatory mechanisms and/or interneuron abnormalities, are detectable even at rest. Given the connectivity strength was relatively low between ipsilateral (LV-LD and RV-RD) and diagonal (LV-RD and RV-LD) region of interest pairs as compared to V-V and D-D (Fig. 2A), caution should be taken in interpreting significant effects involving these pairs (i.e. Table 2, RV-RD/R lesion, LV-RD/D lesion, and Table 3, LD-RV/Below). It may be, however, that synchronous activity between ipsilateral and diagonal grey matter horns occurs in relatively transient intervals at rest, such that the static connectivity over the entire imaging run is low overall. If so, significant effects involving dorsal-ventral horn connectivity may be indicative of relevant physiological processes. Future work using dynamic connectivity analyses, for instance, may shed light on the degree of transient synchrony between ipsilateral and diagonal regions of grey matter in the cord. In terms of the distal effects of lesions on connectivity, a trend was observed in which V-V connectivity was decreased when slices were below a lesion within a descending tract (Fig. 5, P = 0.131). Given ventral grey matter is largely composed of motor neurons that receive input from descending corticospinal axons, it seems reasonable that lesions within these tracts disrupt synchronization of motor neuron populations below the site of pathology. Furthermore, our findings of local and distal effects of structural damage on V-V connectivity are in agreement with the spinal cord injury study by Chen et al. (2015), which showed decreased connectivity within and below the site of injury, but not above. It is important to note, however, that the pathophysiological course of acute spinal cord injury and its subsequent effects on functional connectivity are likely to be qualitatively different from that of demyelinating lesions in multiple sclerosis. The cellular and biochemical cascades accompanied by, for example, gross inflammation and/or breakdown of blood–spinal cord barrier in response to acute injury may have more generalized effects on cord function inferior to injury whereas multiple sclerosis lesion effects, especially during periods of remission, may be more tract-specific (Hausmann, 2003). Nevertheless, the prospects for non-invasive measurement of cord functional connectivity in spinal cord injury and multiple sclerosis patient populations provide an interesting avenue of research for understanding both similarities and differences among dysfunctional cord pathologies. There was a general lack of group level differences in mean functional connectivity metrics between patients with multiple sclerosis and healthy controls (Fig. 2F), especially after controlling for confounding factors (Supplementary Fig. 4). On one hand, it is perhaps encouraging that our connectivity measures are stable across groups and region of interest pairs. However, there are several reasons we may have failed to detect a group difference in mean connectivity values. Our patient sample was relatively heterogeneous in terms of lesion load, clinical disability, and disease duration. This could have potentially conflated early disease stage compensatory increases in connectivity with later stage decreases across the patient group, such that any differences compared to controls were washed out. It is also possible that disease effects could not be teased apart from age effects in our statistical models, given our patient sample was older than controls. Another possibility is that intrinsic spinal cord activity is affected at more acute time windows in relapsing-remitting multiple sclerosis, such as during a relapse phase, and recovers to more normal levels in remitting phases. This would be in-line with both Droby et al. (2016) and Chen et al. (2015), who observed normalization of functional connectivity within several months after acute lesion appearance in the brain and induced injury in cord, respectively. In contrast, the observations from functional MRI activation studies in multiple sclerosis seem somewhat counter to the present findings, where group-level overactivation effects were observed even in relapsing-remitting cohorts. It is important to note that activation in response to a stimulus is not necessarily related to functional connectivity strength, which is instead a measure of the relative synchrony of activity between two regions over time. While functional connectivity patterns at rest are largely analogous to those observed in overt task-based activation studies, resting state networks in the brain are found to be more bilateral in nature. Our results demonstrating distinct, bilateral connectivity between ventral grey matter and between dorsal grey matter in the cord corroborate the general finding of symmetrical networks in the brain. A final caveat to be considered in interpretation of this work as compared to the aforementioned studies of cord activation in multiple sclerosis is that all of these previous studies were performed on 1.5 T scanners and used spin echo functional MRI acquisitions, which are much less sensitive to traditional T2*-weighted BOLD effects and instead are suggested to measure a related, but alternative phenomenon termed SEEP (signal enhancement by extravascular protons) (Stroman et al., 2002, 2003). While the benefits of SEEP functional MRI in the spinal cord have been described (Stroman et al., 2005, 2014), this contrast mechanism still remains controversial among the wider functional MRI community in its ability to measure effects specific to neuronal activation (Jochimsen et al., 2005; Bouwman et al., 2008; Summers and Brooks, 2014). More recent work looking at stimulus-evoked activity and resting state fluctuations in cord, including the present study, have largely transitioned to the employment of BOLD-sensitive acquisitions. Going forward, it will be important to confirm the findings of multiple sclerosis-related overactivation in response to stimulation using BOLD functional MRI. Limitations The present study provides a novel assessment of spinal cord functional connectivity in patients with multiple sclerosis compared to healthy controls, including observations of differential effects of column-specific lesion location on synchronous intrinsic BOLD fluctuations in cervical cord grey matter, and there are a number of limitations to consider in interpreting this work (for a more thorough discussion of some of these limitations, please see the Supplementary material): (i) our rating strategy only provided a coarse identification of lesion location, at the level of columns. Damage to tracts, which are known to specifically innervate cervical cord grey matter (e.g. lateral corticospinal tract), may have more robust effects on connectivity of these regions; (ii) while we took a region of interest-based approach dividing spinal cord grey matter into four horns, connectivity in the cord may be divisible into smaller subnetworks and could benefit from analyses using smaller regions of interest or spatial independent component analysis, which looks for patterns at the voxel level (Kong et al., 2014); (iii) grey matter lesions exist in the spinal cord in multiple sclerosis (Gilmore et al., 2009), as has been appreciated in the brain (Calabrese et al., 2013; Harrison et al., 2015). Future studies should adopt imaging strategies to distinguish healthy and lesioned grey matter in the spinal cord, as this damage may influence connectivity measures (Nelson et al., 2008; De Graaf et al., 2012; Sinnecker et al., 2012); (iv) an important limitation is that, while we looked at the effects of relatively local cervical lesions, structural damage located outside of our imaged volume could have influenced connectivity of the cervical regions analysed in this study. Brain lesions affecting corticospinal tracts, for instance, may modulate spontaneous ventral horn activity or connectivity. With the cervical spine RF-coil used in this study, we were not able to acquire whole-brain images to assess brain lesion burden or location. Future work is needed to understand how brain lesions influence spinal cord function in multiple sclerosis, and a more comprehensive assessment of CNS abnormalities (e.g. inclusion of whole-brain imaging) will be necessary to distinguish upstream/downstream contributions to spinal cord connectivity patterns; (v) it is possible that muscle relaxants and/or sedative drugs such as baclofen or benzodiazepine, which are commonly used to treat patients with spinal cord disorders, have an effect on spinal cord network connectivity. In our patient sample, 7 of 22 subjects were prescribed one of these medications (Supplementary Table 1). We assessed the effect of these medications on mean connectivity and found no significant associations with any region of interest pairings (Supplementary Fig. 5). However, future work with a larger sample size is necessary to confirm this observation and researchers may consider use of these medications as exclusion criteria; and (vi) a final limitation concerns the current practicality of acquiring spinal functional connectivity measurements in a clinical setting, given the relatively complex image processing pipeline used here, which requires some expertise as well as use of a 7 T MRI machine, which are still not widely available. The spinal cord MRI community continues to grow, and along with it is continued development of more user-friendly processing software, similar to what is available for brain imaging. The Spinal Cord Toolbox provides the first comprehensive software package for the processing of spinal cord MRI data, including functional MRI analyses (De Leener et al., 2017). Strategies for robust automatic segmentation of the spinal cord have recently been described as well, which reduce processing time and provide further promise for future clinical application (Asman et al., 2014; De Leener et al., 2016; Prados et al., 2016, 2017; Xu et al., 2016). Regarding 7 T MRI, 510(k) approval for the Siemens 7 T Magnetom Terra system was recently announced by the US Food and Drug Administration, and it is anticipated that this announcement will, in time, facilitate more widespread adoption of 7 T MRI. Nevertheless, several groups have now shown spinal cord functional connectivity is observable at 3 T (e.g. Eippert et al., 2017b; Harita and Stroman, 2017; Weber et al., 2018), suggesting these measures can feasibly be acquired on current clinical scanners. Acknowledgement The authors thank Avyay Reddy for performing quality assurance on preprocessed data. Funding This research was supported by the following funding sources: Novartis IIRP-1456/VUMC 41030 (Pawate), NIH/NIBIB 4T32EB014841-04 (Gore), NIH/NINDS RO1NS092961 (Gore), NIH/NINDS R01NS104149 (Gore), DOD W81XWH-13-1-0073 (Smith), National MS Society, and NIH/NIBIB K99/R00 EB016689 (Barry). The project was also supported by the National Center for Research Resources, grant UL1RR024975-01, which is now at the National Center for Advancing Translational Sciences, grant 2UL1TR000445-06. Supplementary material Supplementary material is available at Brain online. Abbreviation Abbreviation BOLD blood oxygenation level-dependent References Agosta F, Valsasina P, Absinta M, Sala S, Caputo D, Filippi M. Primary progressive multiple sclerosis: tactile-associated functional MR activity in the cervical spinal cord. Radiology  2009; 253: 209– 15. Google Scholar CrossRef Search ADS PubMed  Agosta F, Valsasina P, Caputo D, Stroman PW, Filippi M. Tactile-associated recruitment of the cervical cord is altered in patients with multiple sclerosis. Neuroimage  2008a; 39: 1542– 8. Google Scholar CrossRef Search ADS   Agosta F, Valsasina P, Rocca MA, Caputo D, Sala S, Judica E, et al.   Evidence for enhanced functional activity of cervical cord in relapsing multiple sclerosis. Magn Reson Med  2008b; 59: 1035– 42. Google Scholar CrossRef Search ADS   Asman AJ, Bryan FW, Smith SA, Reich DS, Landman BA. Groupwise multi-atlas segmentation of the spinal cord’s internal structure. Med Image Anal  2014; 18: 460– 71. Google Scholar CrossRef Search ADS PubMed  Barry RL, Gore JC, Smith SA. Development of BOLD-sensitive protocols for imaging the human spinal cord at 7 Tesla. Proc Intl Soc Mag Reson Med  2013; 21: 2278. Barry RL, Rogers BP, Conrad BN, Smith SA, Gore JC. Reproducibility of resting state spinal cord networks in healthy volunteers at 7 Tesla. Neuroimage  2016; 133: 31– 40. Google Scholar CrossRef Search ADS PubMed  Barry RL, Smith SA, Dula AN, Gore JC. Resting state functional connectivity in the human spinal cord. Elife  2014; 3: e02812. Google Scholar PubMed  Barry RL, Strother SC, Gatenby JC, Gore JC. Data-driven optimization and evaluation of 2D EPI and 3D PRESTO for BOLD fMRI at 7 Tesla: I. Focal coverage. Neuroimage  2011; 55: 1034– 43. Google Scholar CrossRef Search ADS PubMed  Barry RL, Vannesjo SJ, By S, Gore JC, Smith SA. Spinal cord MRI at 7T. Neuroimage  2018; 168: 437– 51. Google Scholar CrossRef Search ADS PubMed  Bodini B, Battaglini M, De Stefano N, Khaleeli Z, Barkhof F, Chard D, et al.   T2 lesion location really matters: a 10 year follow-up study in primary progressive multiple sclerosis. J Neurol Neurosurg Psychiatry  2011; 82: 72– 7. Google Scholar CrossRef Search ADS PubMed  Bouwman CJC, Wilmink JT, Mess WH, Backes WH. Spinal cord functional MRI at 3 T: gradient echo echo-planar imaging versus turbo spin echo. Neuroimage  2008; 43: 288– 96. Google Scholar CrossRef Search ADS PubMed  Brown TG. The intrinsic factors in the act of progression in the mammal. Proc R Soc London Ser B  1911; 84: 308– 19. Google Scholar CrossRef Search ADS   Cadotte DW, Stroman PW, Mikulis D, Fehlings MG. A systematic review of spinal fMRI research: outlining the elements of experimental design. J Neurosurgery Spine  2012; 17: 102– 18. Google Scholar CrossRef Search ADS   Calabrese M, Favaretto A, Martini V, Gallo P. Grey matter lesions in MS: from histology to clinical implications. Prion  2013; 7: 20– 7. Google Scholar CrossRef Search ADS PubMed  Chen LM, Mishra A, Yang PF, Wang F, Gore JC. Injury alters intrinsic functional connectivity within the primate spinal cord. Proc Natl Acad Sci USA  2015; 112: 5991– 6. Google Scholar CrossRef Search ADS PubMed  Cohen-Adad J, Wheeler-Kingshott C. Quantitative MRI of the Spinal Cord . New York: Academic Press; 2014. Conrad BN, Maki S, Watchmaker JM, Bailey BA, Barry RL, Smith SA, et al.   BOLD Signal changes in spinal cord with hypercapnia. Proc Intl Soc Mag Reson Med  2017; 25: 5319. Cox RW. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res  1996; 29: 162– 73. Google Scholar CrossRef Search ADS PubMed  De Leener B, Fonov VS, Collins DL, Callot V, Stikov N, Cohen-Adad J. PAM50: unbiased multimodal template of the brainstem and spinal cord aligned with the ICBM152 space. Neuroimage  2018; 165: 170– 9. Google Scholar CrossRef Search ADS PubMed  De Leener B, Lévy S, Dupont SM, Fonov VS, Stikov N, Louis Collins D, et al.   SCT: Spinal Cord Toolbox, an open-source software for processing spinal cord MRI data. Neuroimage  2017; 145: 24– 43. Google Scholar CrossRef Search ADS PubMed  De Leener B, Taso M, Cohen-Adad J, Callot V. Segmentation of the human spinal cord. MAGMA  2016; 29: 125– 53. Google Scholar CrossRef Search ADS PubMed  Droby A, Yuen KSL, Muthuraman M, Reitz SC, Fleischer V, Klein J, et al.   Changes in brain functional connectivity patterns are driven by an individual lesion in MS: a resting-state fMRI study. Brain Imaging Behav  2016; 10: 1117– 26. Google Scholar CrossRef Search ADS PubMed  Dula AN, Pawate S, Dortch RD, Barry RL, George-Durrett KM, Lyttle BD, et al.   Magnetic resonance imaging of the cervical spinal cord in multiple sclerosis at 7T. Mult Scler  2016; 22: 320– 8. Google Scholar CrossRef Search ADS PubMed  Eippert F, Kong Y, Jenkinson M, Tracey I, Brooks JCW. Denoising spinal cord fMRI data: approaches to acquisition and analysis. Neuroimage  2017a; 154: 255– 66. Google Scholar CrossRef Search ADS   Eippert F, Kong Y, Winkler AM, Andersson JL, Finsterbusch J, Büchel C, et al.   Investigating resting-state functional connectivity in the cervical spinal cord at 3T. Neuroimage  2017b; 147: 589– 601. Google Scholar CrossRef Search ADS   Faivre A, Rico A, Zaaraoui W, Crespy L, Reuter F, Wybrecht D, et al.   Assessing brain connectivity at rest is clinically relevant in early multiple sclerosis. Mult Scler J  2012; 18: 1251– 8. Google Scholar CrossRef Search ADS   Filippi M, Rocca MA. Present and future of fMRI in multiple sclerosis. Expert Rev Neurother  2013; 13: 27– 31. Google Scholar CrossRef Search ADS PubMed  Fonov VS, Le Troter A, Taso M, De Leener B, Lévêque G, Benhamou M, et al.   Framework for integrated MRI average of the spinal cord white and gray matter: the MNI-Poly-AMU template. Neuroimage  2014; 102: 817– 27. Google Scholar CrossRef Search ADS PubMed  Gilmore CP, Deluca GC, Bö L, Owens T, Lowe J, Esiri MM, et al.   Spinal cord neuronal pathology in multiple sclerosis. Brain Pathol  2009; 19: 642– 9. Google Scholar CrossRef Search ADS PubMed  Giove F, Garreffa G, Giulietti G, Mangia S, Colonnese C, Maraviglia B. Issues about the fMRI of the human spinal cord. Magn Reson Imaging  2004; 22: 1505– 16. Google Scholar CrossRef Search ADS PubMed  Glover GH, Li TQ, Ress D. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Reson Med  2000; 44: 162– 7. Google Scholar CrossRef Search ADS PubMed  De Graaf WL, Zwanenburg JJM, Visser F, Wattjes MP, Pouwels PJW, Geurts JJG, et al.   Lesion detection at seven Tesla in multiple sclerosis using magnetisation prepared 3D-FLAIR and 3D-DIR. Eur Radiol  2012; 22: 221– 31. Google Scholar CrossRef Search ADS PubMed  Harita S, Stroman PW. Confirmation of resting-state BOLD fluctuations in the human brainstem and spinal cord after identification and removal of physiological noise. Magn Reson Med  2017; 78: 2149– 56. Google Scholar CrossRef Search ADS PubMed  Harrison DM, Roy S, Oh J, Izbudak I, Pham D, Courtney S, et al.   Association of cortical lesion burden on 7-T magnetic resonance imaging with cognition and disability in multiple sclerosis. JAMA Neurol  2015; 72: 1004– 12. Google Scholar CrossRef Search ADS PubMed  Hausmann ON. Post-traumatic inflammation following spinal cord injury. Spinal Cord  2003; 41: 369– 78. Google Scholar CrossRef Search ADS PubMed  Jochimsen TH, Norris DG, Möller HE. Is there a change in water proton density associated with functional magnetic resonance imaging? Magn Reson Med  2005; 53: 470– 3. Google Scholar CrossRef Search ADS PubMed  Kandel ER, Schwartz JH, Jessell TM, Siegelbaum SA, Hudspeth AJ, Mack S. Principles of neural science . 5th edn. New York: McGraw-Hill Companies, Incorporated; 2013. Katura T, Tanaka N, Obata A, Sato H, Maki A. Quantitative evaluation of interrelations between spontaneous low-frequency oscillations in cerebral hemodynamics and systemic cardiovascular dynamics. Neuroimage  2006; 31: 1592– 600. Google Scholar CrossRef Search ADS PubMed  Kearney H, Miller DH, Ciccarelli O. Spinal cord MRI in multiple sclerosis-diagnostic, prognostic and clinical value. Nat Rev Neurol  2015; 11: 327– 38. Google Scholar CrossRef Search ADS PubMed  Kiviniemi V, Ruohonen J, Tervonen O. Separation of physiological very low frequency fluctuation from aliasing by switched sampling interval fMRI scans. Magn Reson Imaging  2005; 23: 41– 6. Google Scholar CrossRef Search ADS PubMed  Kong Y, Eippert F, Beckmann CF, Andersson J, Finsterbusch J, Büchel C, et al.   Intrinsically organized resting state networks in the human spinal cord. Proc Natl Acad Sci USA  2014; 111: 18067– 72. Google Scholar CrossRef Search ADS PubMed  Liu X, Zhou F, Li X, Qian W, Cui J, Zhou IY, et al.   Organization of the intrinsic functional network in the cervical spinal cord: a resting state functional MRI study. Neuroscience  2016; 336: 30– 8. Google Scholar CrossRef Search ADS PubMed  Maieron M, Iannetti GD, Bodurka J, Tracey I, Bandettini PA, Porro CA. Functional responses in the human spinal cord during willed motor actions: evidence for side- and rate-dependent activity. J Neurosci  2007; 27: 4182– 90. Google Scholar CrossRef Search ADS PubMed  Mateo C, Knutsen PM, Tsai PS, Shih AY, Kleinfeld D. Entrainment of arteriole vasomotor fluctuations by neural activity is a basis of blood-oxygenation-level-dependent ‘resting-state’ connectivity. Neuron  2017; 96: 936– 48.e3. Google Scholar CrossRef Search ADS PubMed  McCrea DA, Rybak IA. Organization of mammalian locomotor rhythm and pattern generation. Brain Res Rev  2008; 57: 134– 46. Google Scholar CrossRef Search ADS PubMed  Nelson F, Poonawalla A, Hou P, Wolinsky JS, Narayana PA. 3D MPRAGE improves classification of cortical lesions in multiple sclerosis. Mult Scler  2008; 14: 1214– 19. Google Scholar CrossRef Search ADS PubMed  Prados F, Ashburner J, Blaiotta C, Brosch T, Carballido-Gamio J, Cardoso MJ, et al.   Spinal cord grey matter segmentation challenge. Neuroimage  2017; 152: 312– 29. Google Scholar CrossRef Search ADS PubMed  Prados F, Cardoso MJ, Yiannakas MC, Hoy LR, Tebaldi E, Kearney H, et al.   Fully automated grey and white matter spinal cord segmentation. Sci Rep  2016; 6: 36151. Google Scholar CrossRef Search ADS PubMed  Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: Sensitivity encoding for fast MRI. Magn Reson Med  1999; 42: 952– 62. Google Scholar CrossRef Search ADS PubMed  Purves D, Augustine G, Fitzpatrick D, Hall W, Lamantia AS, White L. Neuroscience . 5th edn. Sunderland, MA: Sinauer Associates, Incorporated; 2012. Rocca MA, Absinta M, Valsasina P, Copetti M, Caputo D, Comi G, et al.   Abnormal cervical cord function contributes to fatigue in multiple sclerosis. Mult Scler J  2012a; 18: 1552– 9. Google Scholar CrossRef Search ADS   Rocca MA, Amato MP, De Stefano N, Enzinger C, Geurts JJ, Penner IK, et al.   Clinical and imaging assessment of cognitive dysfunction in multiple sclerosis. Lancet Neurol  2015; 14: 302– 417. Google Scholar CrossRef Search ADS PubMed  Rocca MA, Valsasina P, Martinelli V, Misci P, Falini A, Comi G, et al.   Large-scale neuronal network dysfunction in relapsing-remitting multiple sclerosis. Neurology  2012b; 79: 1449– 57. Google Scholar CrossRef Search ADS   Rogers BP, Gore JC. Empirical comparison of sources of variation for FMRI connectivity analysis. PLoS One  2008; 3: e3708. Google Scholar CrossRef Search ADS PubMed  Roosendaal SD, Schoonheim MM, Hulst HE, Sanz-Arigita EJ, Smith SM, Geurts JJG, et al.   Resting state networks change in clinically isolated syndrome. Brain  2010; 133: 1612– 21. Google Scholar CrossRef Search ADS PubMed  Sacco R, Bonavita S, Esposito F, Tedeschi G, Gallo A. The contribution of resting state networks to the study of cortical reorganization in MS. Mult Scler Int  2013; 2013: 857807. Google Scholar PubMed  San Emeterio Nateras O, Yu F, Muir ER, Bazan C, Franklin CG, Li W, et al.   Intrinsic resting-state functional connectivity in the human spinal cord at 3.0 T. Radiology  2016; 279: 262– 8. Google Scholar CrossRef Search ADS PubMed  Sinnecker T, Dörr J, Pfueller CF, Harms L, Ruprecht K, Jarius S, et al.   Distinct lesion morphology at 7-T MRI differentiates neuromyelitis optica from multiple sclerosis. Neurology  2012; 79: 708– 14. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Bosma RL, Tsyben A. Somatotopic arrangement of thermal sensory regions in the healthy human spinal cord determined by means of spinal cord functional MRI. Magn Reson Med  2012; 68: 923– 31. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Kornelsen J, Lawrence J, Malisza KL. Functional magnetic resonance imaging based on SEEP contrast: response function and anatomical specificity. Magn Reson Imaging  2005; 23: 843– 50. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Krause V, Malisza KL, Frankenstein UN, Tomanek B. Extravascular proton-density changes as a non-BOLD component of contrast in fMRI of the human spinal cord. Magn Reson Med  2002; 48: 122– 7. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Tomanek B, Krause V, Frankenstein UN, Malisza KL. Functional magnetic resonance imaging of the human brain based on signal enhancement by extravascular protons (SEEP fMRI). Magn Reson Med  2003; 49: 433– 9. Google Scholar CrossRef Search ADS PubMed  Stroman PW, Wheeler-Kingshott C, Bacon M, Schwab JM, Bosma R, Brooks J, et al.   The current state-of-the-art of spinal cord imaging: methods. Neuroimage  2014; 84: 1070– 81. Google Scholar CrossRef Search ADS PubMed  Summers PE, Brooks JCW. Chapter 4.1—Spinal cord fMRI. In: Quantitative MRI of the spinal cord . New York: Academic Press; 2014. p. 221– 39. Tong Y, Lindsey KP, deB Frederick B. Partitioning of physiological noise signals in the brain with concurrent near-infrared spectroscopy and fMRI. J Cereb Blood Flow Metab  2011; 31: 2352– 62. Google Scholar CrossRef Search ADS PubMed  Valsasina P, Agosta F, Absinta M, Sala S, Caputo D, Filippi M. Cervical cord functional MRI changes in relapse-onset MS patients. J Neurol Neurosurg Psychiatry  2010; 81: 405– 8. Google Scholar CrossRef Search ADS PubMed  Valsasina P, Rocca MA, Absinta M, Agosta F, Caputo D, Comi G, et al.   Cervical cord FMRI abnormalities differ between the progressive forms of multiple sclerosis. Hum Brain Mapp  2012; 33: 2072– 80. Google Scholar CrossRef Search ADS PubMed  Weber KA, Sentis AI, Bernadel-Huey ON, Chen Y, Wang X, Parrish TB, et al.   Thermal stimulation alters cervical spinal cord functional connectivity in humans. Neuroscience  2018; 369: 40– 50. Google Scholar CrossRef Search ADS PubMed  Wheeler-Kingshott CA, Stroman PW, Schwab JM, Bacon M, Bosma R, Brooks J, et al.   The current state-of-the-art of spinal cord imaging: applications. Neuroimage  2014; 84: 1082– 93. Google Scholar CrossRef Search ADS PubMed  Winder AT, Echagarruga C, Zhang Q, Drew PJ. Weak correlations between hemodynamic signals and ongoing neural activity during the resting state. Nat Neurosci  2017; 20: 1761– 9. Google Scholar CrossRef Search ADS PubMed  Wu TL, Wang F, Mishra A, Wilson GHIII, Byun N, Chen LM, et al.   Resting-state functional connectivity in the rat cervical spinal cord at 9.4 T. Magn Reson Med  2018; 79: 2773– 83. Google Scholar CrossRef Search ADS PubMed  Xu Z, Conrad BN, Baucom RB, Smith SA, Poulose BK, Landman BA. Abdomen and spinal cord segmentation with augmented active shape models. J Med Imaging  2016; 3: 36002. Google Scholar CrossRef Search ADS   Zhou F, Zhuang Y, Gong H, Wang B, Wang X, Chen Q, et al.   Altered inter-subregion connectivity of the default mode network in relapsing remitting multiple sclerosis: a functional and structural connectivity study. PLoS One  2014; 9: e101198. Google Scholar CrossRef Search ADS PubMed  © The Author(s) (2018). Published by Oxford University Press on behalf of the Guarantors of Brain. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

BrainOxford University Press

Published: Apr 10, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off