- Email: [email protected]
Advances in functional magnetic resonance imaging of the human brainstem
NeuroImage 86 (2014) 91–98
Contents lists available at ScienceDirect
NeuroImage journal homepage: www.elsevier.com/locate/ynimg
Advances in functional magnetic resonance imaging of the human brainstem☆ Florian Beissner a,b,⁎, Andy Schumann a, Franziska Brunn a, Daniela Eisenträger a, Karl-Jürgen Bär a a b
Pain & Autonomics — Integrative Research (PAIR), Department of Psychiatry and Psychotherapy, University Hospital Jena, Germany Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Charlestown, MA, USA
a r t i c l e
i n f o
Article history: Accepted 30 July 2013 Available online 9 August 2013 Keywords: Brainstem Functional magnetic resonance imaging Brainstem nuclei Resting state Physiological noise correction Independent component analysis Intrinsic connectivity networks
a b s t r a c t The brainstem is of tremendous importance for our daily survival, and yet the functional relationships between various nuclei, their projection targets, and afferent regulatory areas remain poorly characterized. The main reason for this lies in the sub-optimal performance of standard neuroimaging methods in this area. In particular, fMRI signals are much harder to detect in the brainstem region compared to cortical areas. Here we describe and validate a new approach to measure activation of brainstem nuclei in humans using standard fMRI sequences and widely available tools for statistical image processing. By spatially restricting an independent component analysis to an anatomically defined brainstem mask, we excluded those areas from the analysis that were strongly affected by physiological noise. This allowed us to identify for the first time intrinsic connectivity networks in the human brainstem and to map brainstem–cortical connectivity purely based on functionally defined regions of interest. © 2013 The Authors. Published by Elsevier Inc. All rights reserved.
Introduction The brainstem is one of the most complicated anatomical entities of the human body. Despite its tremendous importance for our daily survival, the functional relationships between various brainstem nuclei, their projection targets, and afferent regulatory areas remain poorly characterized. The main reason for this lies in the sub-optimal performance of standard neuroimaging methods in this region, including functional magnetic resonance imaging (fMRI) and positron emission tomography (PET). While the application of PET is limited by its spatial resolution, fMRI mainly suffers from an elevated level of physiological noise encountered in the brainstem (Beissner et al., 2011). This noise stems from pulsatile motion of large arteries in the direct vicinity of the brainstem as well as from the flow of cerebro-spinal fluid (CSF) (Harvey et al., 2008; Klose et al., 2000). Both cause strong motion artifacts, especially in the lower brainstem, which cannot be corrected by standard methods, such as realignment. Therefore, BOLD signals in the brainstem region are much harder to detect than those in cortical areas due to a greatly reduced signal-to-noise ratio. While there have been some successful attempts to measure the activity of single brainstem nuclei (D'Ardenne et al., 2008; Eippert et al., 2009; Thompson et al., 2006), the study of inter-nuclear or nucleo-cortical connectivity is still in its infancy. This is due to the fact that all attempts ☆ This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike License, which permits non-commercial use, distribution, and reproduction in any medium, provided the original author and source are credited. ⁎ Corresponding author at: Pain & Autonomics — Integrative Research (PAIR), Department of Psychiatry and Psychotherapy, University Hospital Jena, Germany. E-mail address: fl[email protected] (F. Beissner).
to apply independent component analysis (ICA) (McKeown et al., 1998), one of the standard methods for functional connectivity assessment, to the brainstem, have thus far been unsuccessful. Here we describe and validate a fundamentally different approach to measure activation of brainstem nuclei as well as nucleo-cortical connectivity using a standard fMRI sequence and widely available tools for multivariate statistical image analysis. While most noise-suppression methods, like low-pass filtering and physiological noise regression (Figs. 1d+e), work in the temporal domain, our approach relies mainly on spatial characteristics of physiological noise. The fact that the major part of the noise stems from an area directly adjacent to the brainstem but not from the brainstem itself allows one to exclude those areas from the analysis that are subject to this influence. This is done by restricting the analysis to an anatomically-defined mask of the brainstem. A similar approach has been successfully applied to the cortex before (Formisano et al., 2004). Noise suppression by this masked ICA (“mICA”) approach was followed by source localization by means of ICA, which can be used to detect intrinsic and extrinsic functional connectivity of the brainstem. Another reason to use a brainstem mask is that it prevents results from being driven by the much stronger signals of surrounding subcortical and cerebellar structures. Masked independent component analysis also offers a straightforward approach to measure functional connectivity between the brainstem and cortical areas, avoiding the usual problem of physiological noise interference. To validate our novel approach, we study resting state connectivity networks in the brainstem and show that the majority of them are highly reproducible. Furthermore, we map brainstem–cortical connectivity to identify the brainstem components as neuronal or noise-related based on their projection targets. Our method significantly advances
1053-8119/$ – see front matter © 2013 The Authors. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.neuroimage.2013.07.081
92
F. Beissner et al. / NeuroImage 86 (2014) 91–98
Fig. 1. Conceptual overview of the masked independent component analysis (mICA) approach and comparison of mICA results (f) to those obtained with state-of-the-art methods (d+e). For all analyses, a cuboid containing the brainstem was cut out of the whole-brain dataset (a). b shows the distribution and main sources of physiological noise in a sagittal slice of this cuboid. The results of applying standard ICA to the data without prior noise suppression (c) after regression of physiological signals (d), and after temporal low-pass filtering (e) are shown in the form of the three independent components with the highest uniquely explained variance. Due to their non-Gaussian structure, ICA mainly detected noise components, as can be seen by the peculiar shape of the activations (c–e) coinciding with areas of high physiological noise (b). In stark contrast, the application of an anatomical brainstem mask to exclude areas of high physiological noise leads to ICA results showing activations of individual brainstem nuclei and nuclear complexes (f).
the current capacity to understand the human brainstem and its interactions with other brain regions. Materials and methods Subjects We recruited a sample of 143 healthy subjects (67 males) consisting of students of the local university as well as staff from the university hospital. 43 subjects had to be excluded (30 due to insufficient quality of ECG or respiratory recordings, 1 due to missing MRI data, and 12 due to excessive motion, i.e. N1 mm peak-to-peak, during the measurement). The remaining 100 subjects (52 males) had a mean age of 25.2 ± 9.0 years (mean ± s.d.) and were of normal weight (BMI: 22.8 ± 2.5). None of the subjects had a history of trauma or any other interfering disease. The study was conducted according to the Declaration of Helsinki and all participants gave written informed consent following the guidelines of the local ethics committee who had approved the study. For reproducibility analysis two sub-samples of n = 50 were formed that were matched by age, sex, handedness, body mass index and relative mean motion during the scan. In the following, we will refer to them as discovery and confirmation sample. fMRI measurements All measurements were taken on a 3 T whole body MR scanner (MAGNETOM Trio Tim, Siemens Medical Solutions, Erlangen, Germany) with a 12-channel head matrix coil. The whole measurement consisted of a resting state scan followed by a structural scan. Subjects were asked to keep their eyes closed during the whole measurement. The sequence used for the functional run was gradient-echo echo-planar imaging (GRE-EPI) accelerated by parallel imaging using GRAPPA (Griswold et al., 2002). The parameters were: TE = 30 ms, TR = 2.52 s, GRAPPA factor = 2, PE direction = anterior–posterior, FOV = 220 × 210 mm2, matrix size = 88 × 84, in-plane resolution =
2.5 × 2.5 mm2, slice thickness = 2.5 mm, inter-slice gap = 0.625 mm, and slice tilt = 40°. 45 slices were acquired in ascending order for whole brain coverage including the lower brainstem. The measurement run consisted of 240 volumes with a total length of 10 min and 5 s. The T1-weighted anatomical scan was an MPRAGE with the following parameters: TE = 3.03 ms, TR = 2.3 s, TI = 900 ms, partition thickness = 1 mm, matrix size = 256 × 256, FOV = 280 × 280 mm2, 192 slices, and in-plane resolution = 1.09 × 1.09 mm2. Physiological recordings Electrocardiogram (ECG) and respiration (RESP) were recorded during the MRI scans using an MR-compatible BIOPAC MP150 polygraph (BIOPAC Systems Inc., Goleta, CA, USA). ECG electrodes were arranged in a modified Einthoven's triangle. The sampling rate was 500 Hz for all channels. To remove MRI-related artifacts, ECG signals were bandpass filtered (cutoff: 0.05–35 Hz). The RESP signal was temporally smoothed over 100 samples followed by a simple detection of local maxima. R-waves were extracted after signal decomposition using Daubechies' wavelet of 14th order. Most rapid changes of the ECG signal were detected by thresholding excess of signal components at the highest decomposition level. Resulting inter-beat time series were post-processed by an adaptive filter algorithm described in detail by Wessel (2000). All results were checked off-line by visual inspection. Data preprocessing Software All data processing was carried out using tools from SPM8 (Wellcome Department of Imaging Neuroscience, UCL, London, UK, available at http://www.fil.ion.ucl.ac.uk/spm), FSL5.0 (Oxford Centre for Functional MRI of the Brain, Oxford, UK, available at http://www. fmrib.ox.ac.uk/fsl/), as well as home-written scripts in MATLAB (MathWorks, Natick, MA, USA).
F. Beissner et al. / NeuroImage 86 (2014) 91–98
Anatomical scans Anatomical images were first segmented into gray matter, white matter and cerebrospinal fluid (CSF) using SPM's segment algorithm (Ashburner and Friston, 2005). Gray and white matter maps were added and used as masks to create a brain-extracted version of the anatomical images. Images were then normalized to the non-linear and non-symmetrized version of the ICBM152 (http://www.bic.mni.mcgill. ca/ServicesAtlases/ICBM152NLin2009) using FSL FLIRT und FNIRT and a resolution of 1 mm isotropic for brainstem and 2 mm isotropic for whole-brain analyses. For the brainstem analyses, the images were further cropped to include only the brainstem and adjacent structures in a cuboid caudally confined by the level of the pyramidal decussation and cranially by the AC/PC line. The ventral border was defined by the anterior commissure, the dorsal by the tectum mesencephali. Functional scans Functional images were first corrected for head motion by realignment of each volume to the first volume of the run using FSL's MCFLIRT. Subjects with peak-to-peak motion exceeding 1 mm in any direction were excluded from further analysis. Images were then brain extracted using FSL BET and high-pass filtered with a 100 s cutoff using FSLMATHS followed by regression of physiological signals (see below). The filtered images were linearly coregistered to the anatomical scan using FLIRT with boundary-based registration (Greve and Fischl, 2009) and non-linearly transformed to the ICBM152 using the transformation estimated from the anatomical scan. For brainstem analyses, images were cropped and upsampled to 1 mm isotropic. Spatial smoothing was omitted for these data, as the upsampling from 2.5 to 1 mm generates a sufficiently smooth version. For whole-brain analyses (i.e. back-projection), cropping was omitted and images were upsampled to 2 mm isotropic resolution followed by spatial smoothing with a 5 mm isotropic kernel. Physiological noise regression In addition to our spatial noise reduction approach, we applied regression of physiological signals using FSL's ‘Physiological noise modeling’, PNM (Brooks et al., 2008). Regressors were created for cardiac and respiratory signals using sine and cosine terms of their principal frequencies θC and θR, as well as their first three harmonics (2θC to 4θC and 2θR to 4θR). In addition, interactions of cardiac and respiratory signals were modeled using multiplicative terms of the form sin or cos (AθC ± BθC) with A,B = 1,2 (Kong et al., 2012). Altogether, 32 nuisance regressors were included and treated as voxel-wise confounds in FSL FEAT. As PNM is tailored for application in a general linear model, we had to choose a regressor of interest for the regression step. As low frequencies are of no importance in fMRI experiments, we chose a linear ramp function for this. Masked independent component analysis (mICA) Pre-processed fMRI data of the discovery sample (n = 50) were temporally concatenated and analyzed by probabilistic independent component analysis, pICA (Beckmann and Smith, 2004) using MELODIC 3.13 (FSL5.0). To suppress physiological noise in the spatial domain, pICA was restricted to the volume inside an anatomical brainstem mask, thus, excluding the most problematic areas adjacent to the brainstem from the analysis (masked ICA, mICA, Fig. 1f). The brainstem mask was defined based on gray and white matter tissue maps supplied with the ICBM152 template. The tissue probability was thresholded at 0.9 and remaining non-brainstem regions were removed manually (Supplementary Fig. 1). The relatively high threshold was chosen to result in a smaller mask in order to remove as many CSF-containing (and therefore noisy) voxels from the mask as possible. The final mask comprised the entire medulla, pons, and mesencephalon. Data were projected into a 37-dimensional subspace using probabilistic principal component analysis after voxel-wise de-meaning of the
93
data and normalization of the voxel-wise variance. The decomposition dimensionality was based on the median of all single subject dimensionality estimates using the Laplace approximation to the Bayesian evidence of the model order (Beckmann and Smith, 2004; Minka, 2000). Note that we did not use dimensionality estimates from the whole sample of n = 50, as it has been shown that estimates from analytic methods are strongly influenced by sample size (a higher number of observations leading to higher dimensionality estimates, although the underlying dimensionality stays the same) (Cordes and Nandy, 2006). The whitened observations were decomposed into time-courses and maps by optimizing for non-Gaussian spatial distributions using a fixedpoint iteration technique (Hyvärinen, 1999). Estimated component maps were divided by the standard deviation of the residual noise and thresholded by fitting a mixture model to the histogram of intensity (Beckmann and Smith, 2004). The analysis was repeated with the data of the confirmation sample (n = 50). Comparing mICA to existing methods To illustrate the performance of mICA compared to standard ICA in combination with other noise suppression methods, four separate ICAs were run on the cuboid containing the brainstem (Fig. 1): (1) Unmasked pICA of data without correction for physiological noise, (2) unmasked pICA after regression of physiological signals using PNM, (3) unmasked pICA after low-pass filtering with a 0.1 Hz cutoff, and (4) mICA after regression of physiological signals using PNM. Dimensionalities were estimated on the group level for each analysis using the Laplace approximation to the Bayesian evidence of the model order. The three independent components with the highest uniquely explained variance in each approach were extracted for visual comparison. Furthermore, variance maps were calculated from each of the four functional datasets. Reproducibility and symmetry analysis of mICA results Reproducibility analysis of the brainstem ICN was accomplished by calculating spatial cross-correlation coefficients for all possible combinations of the un-thresholded independent components from the discovery and the confirmation sample. The resulting cross-correlation matrix was then sorted using the Hungarian sorting algorithm for a linear assignment problem (Munkres, 1957). This yielded a one-to-one mapping maximizing the sum of all coefficients. Components with a correlation coefficient above 0.5 were considered reproducible. Symmetry of the brainstem ICN in the discovery sample was assessed following an approach similar to reproducibility analysis. Correlation coefficients were calculated for all components of the discovery sample and their left–right flipped counterparts. Again, Hungarian sorting yielded a one-to-one mapping, where components with an intrinsic symmetry were mapped on themselves and all other components on their best-matching symmetric partners. Components with a correlation coefficient larger than 0.5 were considered intrinsically symmetric and symmetric partners, respectively. Back-projection analysis to study cortical connectivity and specificity of brainstem components We back-projected the associated time-courses of all 37 mICA components of the discovery sample onto the whole-brain functional data (Calhoun et al., 2001) using a general linear model (Fig. 2). This was done to find all brain areas exhibiting a temporal behavior similar to that of the brainstem components under consideration. The reason for this was twofold. Firstly, we wanted to differentiate specific components (i.e. those of neuronal origin) from unspecific ones (i.e. those related to noise) (Fig. 2). Secondly, we were interested in the cortical and subcortical functional connectivity profile of the brainstem
94
F. Beissner et al. / NeuroImage 86 (2014) 91–98
proportions. Specificity was defined as s = 1-pCSF. Components were considered specific, if their s-value was within one standard deviation of the mean s. mICA vs. mICA + PNM To assess the impact of physiological noise regression (PNM) as a pre-processing step, when using the masked ICA approach, we also repeated mICA without prior noise regression for both discovery and confirmation sample. For the sake of comparability, the number of components was kept at 37. Matching components between the two samples (with and without PNM) were once more identified by calculating spatial cross-correlation coefficients and subsequent Hungarian sorting. We assessed reproducibility of components derived without prior PNM using the same approach as described for the PNM sample. A two-sided two-sample t-test was calculated for the hypothesis of reproducibility differences between the two analyses. A p-value of 0.05 was considered significant. Results mICA compared to existing methods
Fig. 2. Conceptual overview of the approach used to assess specificity of the brainstem components (see Fig. 4 for results). For each component, the associated time-course was back-projected into the whole-brain data to identify regions exhibiting a similar timecourse. The sum of these regions can be interpreted as the cortical (and sub-cortical) connectivity profile of the brainstem component. In a second step the thresholded connectivity patterns were scalarly multiplied with a probability mask of cerebro-spinal fluid (CSF) to calculate the percentage of activated voxels that lay within CSF regions (pCSF). Specificity was defined as s = 1-pCSF.
components. The analysis was done separately for each subject by splitting the concatenated time-courses of the ICA into single-subject timecourses. The single-subject parameter estimates were then tested for statistically significant areas using a non-parametric permutation test implemented in RANDOMISE 2.9 (FSL5.0). 5000 permutations of the data were calculated. A voxel-wise p-value of 0.01, corrected for multiple comparisons by controlling the family-wise error (FWE), was considered significant. For specificity analysis, we calculated the percentage pCSF of the overall volume of activated voxels that fell into cerebro-spinal fluid (CSF) regions as compared to white and gray matter. The rationale for this was that a high value would indicate a strong contribution of CSF signal to the brainstem component, supporting a non-neuronal origin. The thresholded, binarized p-maps were scalarly multiplied with the tissue probability map of CSF distribution available for the ICBM152 template. The same procedure was repeated for white and gray matter maps. Prior to multiplication all probability maps had been thresholded at a tissue probability of 0.5. pCSF was defined as the CSF proportion of activated voxels normalized by the sum of the white and gray matter
The results of the analysis under the three different conditions (without physiological noise correction, after regression of physiological signals, and after temporal low-pass filtering) are shown in Fig. 1. Dimensionality estimates were 16 for no correction, 22 for physiological noise regression, and 55 for low-pass filtering. In contrast, the dimensionality estimate for the mICA approach was 95. For each analysis, the three independent components with the highest uniquely explained variance are shown, while the entire set of components can be found in Supplementary Figs. 2, 3, and 4 (for the alternative approaches). The dominant independent components of all three alternative analyses showed similar artifacts, which could easily be sub-divided into two classes. The majority of ICA components showed a pattern of activity that was chiefly located outside the brainstem. A comparison with physiological noise patterns (Fig. 1b) revealed a strong similarity between their spatial distribution and that of the ICA results. This is expected, as ICA algorithms identify signal sources by their nonGaussianity, a property shared by almost all physiological noise signals. As noise signals have stronger signal intensities than BOLD signals, they dominate the ICA results. The second class of artifacts was also found in the results of all three analyses: Extensive activation clusters that were partly located in the brainstem, but primarily comprised subcortical and cerebellar regions. These components result from BOLD signals in nonbrainstem regions whose signal intensities exceed those within the brainstem. Although they are valid components depicting neuronal activity, they cannot be used to study the behavior of brainstem nuclei, as they usually form large clusters (cf. Supplementary Fig. 3, lower left corner). When applying an anatomical brainstem mask, and restricting the ICA within that mask (i.e. the mICA approach), the analysis yielded starkly contrasting results: as shown in Fig. 1f, activation clusters of single nuclei and nuclear complexes were made evident, and a relative lack of artifacts could be observed (Fig. 1f). Reproducibility and symmetry of mICA results The results of the mICA analysis for the discovery sample are shown in Fig. 3. For comparison, those of the confirmation sample can be found in Supplementary Fig. 5. Components were ordered by their anatomical localization from rostral to caudal and show a rather equal distribution between mesencephalon (10 ICs), pons (20 ICs), and medulla (7 ICs) taking into account the relative volume of these parts of the brainstem.
F. Beissner et al. / NeuroImage 86 (2014) 91–98
95
Fig. 3. Results of a 37-dimensional masked independent component analysis of 50 healthy subjects. Brainstem components are presented at a mixture model threshold of 0.5 and ordered in rostro-caudal direction. Reproducibility (rep) was assessed by spatial cross-correlation with components in a confirmation sample (Supplementary Fig. 5), while specificity (spec) scores are based on whole-brain connectivity of the components (Fig. 2). Graphs on the lower right show the distributions of these scores as well as the cutoff used to identify non-reproducible (ο) and non-specific (i.e. noise-related) components (#). Note that the ordering of the components is different in both graphs.
The majority (32/37) of all components in the discovery sample were reproducible, i.e. they had a correlation coefficient higher than r = 0.5 with their best matching partner from the confirmation sample. The average correlation coefficient of the reproducible components was 0.65 ± 0.11 (mean ± STD). Components showed a high symmetry, in that 32 of the components were either intrinsically symmetric (14 components) or had a contralateral partner (18 components). This means that all these components had a correlation coefficient higher than r = 0.5 with one of the other components of the sample. The average correlation coefficient of the symmetric components was 0.64 ± 0.08 (mean ± STD).
Whole-brain connectivity and specificity of brainstem components Whole-brain connectivity profiles showed a large variability across brainstem components. A group of representative components is shown in Fig. 4. Some of them resembled known large-scale brain networks, like the default mode network (Fig. 4a), while others did not (Figs. 4b+c). Three classes of components could be distinguished by eye: The first contained the abovementioned networks of cortical and subcortical structures, where connectivity was clearly dominated by gray matter regions (Figs. 4a–c). In the second class, profiles comprised a much higher percentage of CSF and white matter regions (Figs. 4d–f).
96
F. Beissner et al. / NeuroImage 86 (2014) 91–98
Fig. 4. Results of the cortical connectivity/specificity analysis (see Fig. 2 for method). Three representative components are shown for the specific and unspecific group together with their specificity values (spec). While the cortical connectivity profiles of the specific components (a–c) are dominated by gray matter regions, unspecific components (d–f) exhibit noticeable connectivity with CSF regions. An interesting observation that further corroborates the validity of our approach is the median pontine component in (a) showing a connectivity profile that largely resembles the default mode network, a connection that has been reported before (Habas et al., 2009).
Specificity analysis confirmed this distinction, as the latter components had higher pCSF and therewith lower specificity values. The third class comprised only two components, namely the caudal and rostral periaqueductal gray, which showed only very localized connectivity that was, however, dominated by CSF of the aqueduct leading to their classification as non-specific. mICA vs. mICA + PNM Omitting physiological noise regression in the mICA analysis yielded results that were very similar to that of the original analysis (Supplementary Fig. 6). The average correlation coefficient between both sets of components was 0.79 ± 0.16 (mean ± STD). Only two components had correlation coefficients below 0.5. One of them was a bilateral activation pattern of the superior colliculi, the other one a putative noise component. Reproducibility did not significantly differ between the PNM and non-PNM analyses (0.59 ± 0.21 vs. 0.61 ± 0.21, p = 0.46) (Supplementary Fig. 6). Discussion In this study, we combined a multivariate data-driven analysis method (ICA) with an anatomical mask and applied it to resting-state fMRI data, to study the activity of brainstem nuclei and intrinsic connectivity networks. To the best of our knowledge, it constitutes the first successful attempt to depict brainstem nuclei as well as their connectivity to cortical areas in humans using ICA. The key distinction between our masked ICA approach and previous methods is the radically different treatment of physiological noise by exploiting mainly spatial instead of temporal characteristics.
Instead of using a two-step procedure of mICA and back-projection to arrive at brainstem–cortex connectivity profiles, one could also think of using a whole-brain ICA approach. There are, however, several reasons, why this common approach usually does not identify brainstem foci associated with the cortical component. One is that whole-brain group analyses require at least a minimal amount of spatial smoothing. This pre-processing step inevitably leads to contamination of the brainstem with physiological noise from the adjacent CSFcontaining regions. Furthermore, the optimal amount of smoothing to depict brainstem nuclei may be different from that commonly used for cortical studies (Beissner et al., 2011). The most important reason, why brainstem- and whole-brain-centered approaches produce different results, is that whole-brain ICAs are largely driven by cortical signals due to the much larger relative volume of the cortex compared to the brainstem. Thus, brainstem nuclei are only detected, if they belong to one of the known large-scale resting state networks. Fig. 4a shows such an example, where back-projection of a median pontine component yielded the default mode network — a result that corroborates previous reports of pontine nuclei being part of the default mode network (Habas et al., 2009). In general, however, the cortical counterpart of brainstem networks may be small compared to large-scale brain networks and may, thus, remain completely undetected by whole-brain ICA. As Figs. 4b+c show, such cortical components may indeed look very different from those commonly reported in the literature. Concerning reproducibility and specificity of our results, we found that brainstem resting-state components show high reproducibility, at least across samples acquired at the same MR scanner. However, our results also show that reproducibility must not be equated with specificity. The fact that the five components with the lowest specificity scores all showed moderate reproducibility, underlines that even components
F. Beissner et al. / NeuroImage 86 (2014) 91–98
that are reproducible may still be of non-neuronal, i.e. artifactual, origin. This result, however, is not unexpected, as physiological noise is highly structured and there is no reason to believe that its spatial profile should differ between subjects. As mICA relies mainly on a spatial approach to remove physiological noise (by excluding highly affected areas from the analysis), we have combined it with PNM, a temporal noise reduction method (Brooks et al., 2008). As the two approaches are complementary, one could think that their combination should yield better results than either of it alone. However, while pre-processing with PNM led to a higher dimensionality estimate compared to the uncorrected analysis (Supplementary Figs. 2 vs. 3), we were not able to find significant differences in reproducibility when a fixed dimensionality was used for the comparison (Supplementary Fig. 6). The decision whether or not an additional temporal noise regression step should be applied to the data has to be made individually for each study. On the one hand, a combination of both techniques (PNM + mICA) may be a more conservative approach, whose potential should be assessed in further studies. On the other hand, every regression of physiological “noise” bears the risk of removing meaningful signal from the data (Iacovella and Hasson, 2011; Khalili-Mahani et al., 2012). This is especially relevant for studies of autonomic nuclei or of those triggering autonomic responses (e.g. nociceptive). One great advantage of omitting PNM is of course the applicability of mICA to existing fMRI data sets that lack physiological recordings. However, the authors should take great care to control for possible residual noise. Another pre-processing step that may be optional is the upsampling of the brainstem data to 1 mm. While this step clearly improves post-hoc comparison of functional and structural data, a resolution of 2 mm was found to produce very similar results later (results not shown). Several limitations should be addressed. Firstly, we used a restingstate experiment, whereas specific tasks may have facilitated the identification of nuclei. However, correspondence of the brain's functional architecture during tasks and resting-state has recently been demonstrated for the cortex (Smith et al., 2009), a concept that can probably be transferred to the brainstem. Furthermore, validation of brainstem nuclei based on task-fMRI would be hampered by the lack of tasks that activate only a small number of nuclei. For example, letting subjects move their eyes to a fixed cue would not only activate the abducens nuclei, but also visual, attention, vestibular and many more systems localized side-by-side in the brainstem. Another problem is that almost all nuclei have multiple functions (Koutcherov et al., 2009) and for many of them, these functions are still unknown. Secondly, by restricting the ICA to the brainstem region, we lost spatial information from the rest of the brain that could otherwise be used to identify imaging artifacts (e.g., ghosting and motion). However, parts of the information lost by cropping of the images are regained by our back-projection approach. As Supplementary Figs. 2–4 show, artifacts in the brainstem region also have specific patterns making their visual identification possible. Our results raise the prospect that the measurement of intrinsic and extrinsic connectivity of the brainstem will become standard in fMRI studies, both in humans and in small animal imaging and complement information from cortical and subcortical regions. Studies of brainstem–cortex connectivity in humans will prove crucial for the understanding of functional systems, like the neuromodulatory (Briand et al., 2007) and nociceptive (Heinricher et al., 2009) systems or the central autonomic network (Benarroch, 1993), as all of them have vital nodes in the brainstem, and extensive connections with the cortex. Furthermore, these systems have critical clinical relevance; notable examples include the neuromodulatory systems in major depression (Belmaker and Agam, 2008) and Parkinson's disease (Obeso et al., 2010), the nociceptive system in chronic pain syndromes (Tracey and Mantyh, 2007), and the autonomic nervous system in hypertension and cardiovascular diseases (Brook and Julius, 2000).
97
Finally, it is the hope of the authors that the brainstem will lose its label of a terra incognita and soon become a region of major interest in the neuroimaging community. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2013.07.081.
Acknowledgments FBe thanks Vitaly Napadow for helpful discussions. FBe was supported by the German Research Foundation Grant BE4677/1-1. FBr and AS were supported by the Horst-Görtz-Foundation.
References Ashburner, J., Friston, K.J., 2005. Unified segmentation. NeuroImage 26 (3), 839–851. http://dx.doi.org/10.1016/j.neuroimage.2005.02.018. Beckmann, C.F., Smith, S.M., 2004. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23 (2), 137–152. Beissner, F., Deichmann, R., Baudrexel, S., 2011. fMRI of the brainstem using dual-echo EPI. NeuroImage 55 (4), 1593–1599. http://dx.doi.org/10.1016/j.neuroimage.2011.01.042. Belmaker, R.H., Agam, G., 2008. Major depressive disorder. N. Engl. J. Med. 358 (1), 55–68. http://dx.doi.org/10.1056/NEJMra073096. Benarroch, E., 1993. The central autonomic network: functional organization, dysfunction, and perspective. Mayo Clin. Proc. Mayo Clin. 68, 988. Briand, L.A., Gritton, H., Howe, W.M., Young, D.A., Sarter, M., 2007. Modulators in concert for cognition: modulator interactions in the prefrontal cortex. Prog. Neurobiol. 83 (2), 69–91. http://dx.doi.org/10.1016/j.pneurobio.2007.06.007. Brook, R.D., Julius, S., 2000. Autonomic imbalance, hypertension, and cardiovascular risk. Am. J. Hypertens. 13 (6 Pt 2), 112S–122S. Brooks, J.C.W., Beckmann, C.F., Miller, K.L., Wise, R.G., Porro, C.A., Tracey, I., Jenkinson, M., 2008. Physiological noise modelling for spinal functional magnetic resonance imaging studies. NeuroImage 39 (2), 680–692. http://dx.doi.org/10.1016/j.neuroimage. 2007.09.018. Calhoun, V.D., Adali, T., Pearlson, G.D., Pekar, J.J., 2001. A method for making group inferences from functional MRI data using independent component analysis. Hum. Brain 151, 140–151. http://dx.doi.org/10.1002/hbm. Cordes, D., Nandy, R., 2006. Estimation of the intrinsic dimensionality of fMRI data. NeuroImage 29 (1), 145–154. D'Ardenne, K., McClure, S.M., Nystrom, L.E., Cohen, J.D., 2008. BOLD responses reflecting dopaminergic signals in the human ventral tegmental area. Science 319 (5867), 1264–1267. http://dx.doi.org/10.1126/science.1150605. Eippert, F., Bingel, U., Schoell, E.D., Yacubian, J., Klinger, R., Lorenz, J., Büchel, C., 2009. Activation of the opioidergic descending pain control system underlies placebo analgesia. Neuron 63 (4), 533–543. http://dx.doi.org/10.1016/j.neuron.2009.07.014. Formisano, E., Esposito, F., Di Salle, F., Goebel, R., 2004. Cortex-based independent component analysis of fMRI time series. Magn. Reson. Imaging 22 (10), 1493–1504. http:// dx.doi.org/10.1016/j.mri.2004.10.020. Greve, D.N., Fischl, B., 2009. Accurate and robust brain image alignment using boundarybased registration. NeuroImage 48 (1), 63–72. http://dx.doi.org/10.1016/j.neuroimage. 2009.06.060. Griswold, M.A., Jakob, P.M., Heidemann, R.M., Nittka, M., Jellus, V., Wang, J., Kiefer, B., et al., 2002. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn. Reson. Med. 47 (6), 1202–1210. http://dx.doi.org/10.1002/mrm.10171. Habas, C., Kamdar, N., Nguyen, D., Prater, K., Beckmann, C.F., Menon, V., Greicius, M.D., 2009. Distinct cerebellar contributions to intrinsic connectivity networks. J. Neurosci. 29 (26), 8586–8594. http://dx.doi.org/10.1523/JNEUROSCI.1868-09.2009. Harvey, A.K., Pattinson, K.T.S., Brooks, J.C.W., Mayhew, S.D., Jenkinson, M., Wise, R.G., 2008. Brainstem functional magnetic resonance imaging: disentangling signal from physiological noise. J. Magn. Reson. Imaging 28 (6), 1337–1344. Heinricher, M.M., Tavares, I., Leith, J.L., Lumb, B.M., 2009. Descending control of nociception: specificity, recruitment and plasticity. Brain Res. Rev. 60 (1), 214–225. http://dx.doi.org/10.1016/j.brainresrev.2008.12.009. Hyvärinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Netw. 10 (3), 626–634. Iacovella, V., Hasson, U., 2011. The relationship between BOLD signal and autonomic nervous system functions: implications for processing of “physiological noise”. Magn. Reson. Imaging 29, 1338–1345. Khalili-Mahani, N., Chang, C., Van Osch, M.J., Veer, I.M., Van Buchem, M.A., Dahan, A., Beckmann, C.F., et al., 2012. The impact of “physiological correction” on functional connectivity analysis of pharmacological resting state fMRI. NeuroImage 65, 499–510. http://dx.doi.org/10.1016/j.neuroimage.2012.09.044. Klose, U., Strik, C., Kiefer, C., Grodd, W., 2000. Detection of a relation between respiration and CSF pulsation with an echoplanar technique. J. Magn. Reson. Imaging:JMRI2 11 (4), 438–444. Kong, Y., Jenkinson, M., Andersson, J., Tracey, I., Brooks, J.C.W., 2012. Assessment of physiological noise modelling methods for functional imaging of the spinal cord. NeuroImage 60 (2), 1538–1549. http://dx.doi.org/10.1016/j.neuroimage.2011.11.077. Koutcherov, Y., Huang, X.F., Halliday, G., Paxinos, G.T., 2009. Organization of human brain stem nuclei. In: Paxinos, G.T., Mai, J.K. (Eds.), The Human Nervous System. Elsevier, London, pp. 267–320.
98
F. Beissner et al. / NeuroImage 86 (2014) 91–98
McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S., Bell, A.J., Sejnowski, T.J., 1998. Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapp. 6 (3), 160–188. Minka, T., 2000. Automatic choice of dimensionality for PCA. Technical Report 514. Munkres, J., 1957. Algorithms for the assignment and transportation problems. J. Soc. Ind. Appl. Math. 5 (1), 32–38. Obeso, J.A., Rodriguez-Oroz, M.C., Goetz, C.G., Marin, C., Kordower, J.H., Rodriguez, M., Hirsch, E.C., et al., 2010. Missing pieces in the Parkinson's disease puzzle. Nat. Med. 16 (6), 653–661. http://dx.doi.org/10.1038/nm.2165. Smith, S.M., Fox, P., Miller, K., 2009. Correspondence of the brain's functional architecture during activation and rest. Proceedings of the National Academy of Sciences of
the United States of America 106 (31), 13040–13045. http://dx.doi.org/10.1073/ pnas.0905267106. Thompson, S.K., Von Kriegstein, K., Deane-Pratt, A., Marquardt, T., Deichmann, R., Griffiths, T.D., McAlpine, D., 2006. Representation of interaural time delay in the human auditory midbrain. Nat. Neurosci. 9 (9), 1096–1098. http://dx.doi.org/10.1038/ nn1755. Tracey, I., Mantyh, P.W., 2007. The cerebral signature for pain perception and its modulation. Neuron 55 (3), 377–391. http://dx.doi.org/10.1016/j.neuron.2007.07. 012. Wessel, N., 2000. Nonlinear analysis of complex phenomena in cardiological data. Herzschrittmacherther. Elektrophysiol. 11, 159–173.
Contents lists available at ScienceDirect
NeuroImage journal homepage: www.elsevier.com/locate/ynimg
Advances in functional magnetic resonance imaging of the human brainstem☆ Florian Beissner a,b,⁎, Andy Schumann a, Franziska Brunn a, Daniela Eisenträger a, Karl-Jürgen Bär a a b
Pain & Autonomics — Integrative Research (PAIR), Department of Psychiatry and Psychotherapy, University Hospital Jena, Germany Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Charlestown, MA, USA
a r t i c l e
i n f o
Article history: Accepted 30 July 2013 Available online 9 August 2013 Keywords: Brainstem Functional magnetic resonance imaging Brainstem nuclei Resting state Physiological noise correction Independent component analysis Intrinsic connectivity networks
a b s t r a c t The brainstem is of tremendous importance for our daily survival, and yet the functional relationships between various nuclei, their projection targets, and afferent regulatory areas remain poorly characterized. The main reason for this lies in the sub-optimal performance of standard neuroimaging methods in this area. In particular, fMRI signals are much harder to detect in the brainstem region compared to cortical areas. Here we describe and validate a new approach to measure activation of brainstem nuclei in humans using standard fMRI sequences and widely available tools for statistical image processing. By spatially restricting an independent component analysis to an anatomically defined brainstem mask, we excluded those areas from the analysis that were strongly affected by physiological noise. This allowed us to identify for the first time intrinsic connectivity networks in the human brainstem and to map brainstem–cortical connectivity purely based on functionally defined regions of interest. © 2013 The Authors. Published by Elsevier Inc. All rights reserved.
Introduction The brainstem is one of the most complicated anatomical entities of the human body. Despite its tremendous importance for our daily survival, the functional relationships between various brainstem nuclei, their projection targets, and afferent regulatory areas remain poorly characterized. The main reason for this lies in the sub-optimal performance of standard neuroimaging methods in this region, including functional magnetic resonance imaging (fMRI) and positron emission tomography (PET). While the application of PET is limited by its spatial resolution, fMRI mainly suffers from an elevated level of physiological noise encountered in the brainstem (Beissner et al., 2011). This noise stems from pulsatile motion of large arteries in the direct vicinity of the brainstem as well as from the flow of cerebro-spinal fluid (CSF) (Harvey et al., 2008; Klose et al., 2000). Both cause strong motion artifacts, especially in the lower brainstem, which cannot be corrected by standard methods, such as realignment. Therefore, BOLD signals in the brainstem region are much harder to detect than those in cortical areas due to a greatly reduced signal-to-noise ratio. While there have been some successful attempts to measure the activity of single brainstem nuclei (D'Ardenne et al., 2008; Eippert et al., 2009; Thompson et al., 2006), the study of inter-nuclear or nucleo-cortical connectivity is still in its infancy. This is due to the fact that all attempts ☆ This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike License, which permits non-commercial use, distribution, and reproduction in any medium, provided the original author and source are credited. ⁎ Corresponding author at: Pain & Autonomics — Integrative Research (PAIR), Department of Psychiatry and Psychotherapy, University Hospital Jena, Germany. E-mail address: fl[email protected] (F. Beissner).
to apply independent component analysis (ICA) (McKeown et al., 1998), one of the standard methods for functional connectivity assessment, to the brainstem, have thus far been unsuccessful. Here we describe and validate a fundamentally different approach to measure activation of brainstem nuclei as well as nucleo-cortical connectivity using a standard fMRI sequence and widely available tools for multivariate statistical image analysis. While most noise-suppression methods, like low-pass filtering and physiological noise regression (Figs. 1d+e), work in the temporal domain, our approach relies mainly on spatial characteristics of physiological noise. The fact that the major part of the noise stems from an area directly adjacent to the brainstem but not from the brainstem itself allows one to exclude those areas from the analysis that are subject to this influence. This is done by restricting the analysis to an anatomically-defined mask of the brainstem. A similar approach has been successfully applied to the cortex before (Formisano et al., 2004). Noise suppression by this masked ICA (“mICA”) approach was followed by source localization by means of ICA, which can be used to detect intrinsic and extrinsic functional connectivity of the brainstem. Another reason to use a brainstem mask is that it prevents results from being driven by the much stronger signals of surrounding subcortical and cerebellar structures. Masked independent component analysis also offers a straightforward approach to measure functional connectivity between the brainstem and cortical areas, avoiding the usual problem of physiological noise interference. To validate our novel approach, we study resting state connectivity networks in the brainstem and show that the majority of them are highly reproducible. Furthermore, we map brainstem–cortical connectivity to identify the brainstem components as neuronal or noise-related based on their projection targets. Our method significantly advances
1053-8119/$ – see front matter © 2013 The Authors. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.neuroimage.2013.07.081
92
F. Beissner et al. / NeuroImage 86 (2014) 91–98
Fig. 1. Conceptual overview of the masked independent component analysis (mICA) approach and comparison of mICA results (f) to those obtained with state-of-the-art methods (d+e). For all analyses, a cuboid containing the brainstem was cut out of the whole-brain dataset (a). b shows the distribution and main sources of physiological noise in a sagittal slice of this cuboid. The results of applying standard ICA to the data without prior noise suppression (c) after regression of physiological signals (d), and after temporal low-pass filtering (e) are shown in the form of the three independent components with the highest uniquely explained variance. Due to their non-Gaussian structure, ICA mainly detected noise components, as can be seen by the peculiar shape of the activations (c–e) coinciding with areas of high physiological noise (b). In stark contrast, the application of an anatomical brainstem mask to exclude areas of high physiological noise leads to ICA results showing activations of individual brainstem nuclei and nuclear complexes (f).
the current capacity to understand the human brainstem and its interactions with other brain regions. Materials and methods Subjects We recruited a sample of 143 healthy subjects (67 males) consisting of students of the local university as well as staff from the university hospital. 43 subjects had to be excluded (30 due to insufficient quality of ECG or respiratory recordings, 1 due to missing MRI data, and 12 due to excessive motion, i.e. N1 mm peak-to-peak, during the measurement). The remaining 100 subjects (52 males) had a mean age of 25.2 ± 9.0 years (mean ± s.d.) and were of normal weight (BMI: 22.8 ± 2.5). None of the subjects had a history of trauma or any other interfering disease. The study was conducted according to the Declaration of Helsinki and all participants gave written informed consent following the guidelines of the local ethics committee who had approved the study. For reproducibility analysis two sub-samples of n = 50 were formed that were matched by age, sex, handedness, body mass index and relative mean motion during the scan. In the following, we will refer to them as discovery and confirmation sample. fMRI measurements All measurements were taken on a 3 T whole body MR scanner (MAGNETOM Trio Tim, Siemens Medical Solutions, Erlangen, Germany) with a 12-channel head matrix coil. The whole measurement consisted of a resting state scan followed by a structural scan. Subjects were asked to keep their eyes closed during the whole measurement. The sequence used for the functional run was gradient-echo echo-planar imaging (GRE-EPI) accelerated by parallel imaging using GRAPPA (Griswold et al., 2002). The parameters were: TE = 30 ms, TR = 2.52 s, GRAPPA factor = 2, PE direction = anterior–posterior, FOV = 220 × 210 mm2, matrix size = 88 × 84, in-plane resolution =
2.5 × 2.5 mm2, slice thickness = 2.5 mm, inter-slice gap = 0.625 mm, and slice tilt = 40°. 45 slices were acquired in ascending order for whole brain coverage including the lower brainstem. The measurement run consisted of 240 volumes with a total length of 10 min and 5 s. The T1-weighted anatomical scan was an MPRAGE with the following parameters: TE = 3.03 ms, TR = 2.3 s, TI = 900 ms, partition thickness = 1 mm, matrix size = 256 × 256, FOV = 280 × 280 mm2, 192 slices, and in-plane resolution = 1.09 × 1.09 mm2. Physiological recordings Electrocardiogram (ECG) and respiration (RESP) were recorded during the MRI scans using an MR-compatible BIOPAC MP150 polygraph (BIOPAC Systems Inc., Goleta, CA, USA). ECG electrodes were arranged in a modified Einthoven's triangle. The sampling rate was 500 Hz for all channels. To remove MRI-related artifacts, ECG signals were bandpass filtered (cutoff: 0.05–35 Hz). The RESP signal was temporally smoothed over 100 samples followed by a simple detection of local maxima. R-waves were extracted after signal decomposition using Daubechies' wavelet of 14th order. Most rapid changes of the ECG signal were detected by thresholding excess of signal components at the highest decomposition level. Resulting inter-beat time series were post-processed by an adaptive filter algorithm described in detail by Wessel (2000). All results were checked off-line by visual inspection. Data preprocessing Software All data processing was carried out using tools from SPM8 (Wellcome Department of Imaging Neuroscience, UCL, London, UK, available at http://www.fil.ion.ucl.ac.uk/spm), FSL5.0 (Oxford Centre for Functional MRI of the Brain, Oxford, UK, available at http://www. fmrib.ox.ac.uk/fsl/), as well as home-written scripts in MATLAB (MathWorks, Natick, MA, USA).
F. Beissner et al. / NeuroImage 86 (2014) 91–98
Anatomical scans Anatomical images were first segmented into gray matter, white matter and cerebrospinal fluid (CSF) using SPM's segment algorithm (Ashburner and Friston, 2005). Gray and white matter maps were added and used as masks to create a brain-extracted version of the anatomical images. Images were then normalized to the non-linear and non-symmetrized version of the ICBM152 (http://www.bic.mni.mcgill. ca/ServicesAtlases/ICBM152NLin2009) using FSL FLIRT und FNIRT and a resolution of 1 mm isotropic for brainstem and 2 mm isotropic for whole-brain analyses. For the brainstem analyses, the images were further cropped to include only the brainstem and adjacent structures in a cuboid caudally confined by the level of the pyramidal decussation and cranially by the AC/PC line. The ventral border was defined by the anterior commissure, the dorsal by the tectum mesencephali. Functional scans Functional images were first corrected for head motion by realignment of each volume to the first volume of the run using FSL's MCFLIRT. Subjects with peak-to-peak motion exceeding 1 mm in any direction were excluded from further analysis. Images were then brain extracted using FSL BET and high-pass filtered with a 100 s cutoff using FSLMATHS followed by regression of physiological signals (see below). The filtered images were linearly coregistered to the anatomical scan using FLIRT with boundary-based registration (Greve and Fischl, 2009) and non-linearly transformed to the ICBM152 using the transformation estimated from the anatomical scan. For brainstem analyses, images were cropped and upsampled to 1 mm isotropic. Spatial smoothing was omitted for these data, as the upsampling from 2.5 to 1 mm generates a sufficiently smooth version. For whole-brain analyses (i.e. back-projection), cropping was omitted and images were upsampled to 2 mm isotropic resolution followed by spatial smoothing with a 5 mm isotropic kernel. Physiological noise regression In addition to our spatial noise reduction approach, we applied regression of physiological signals using FSL's ‘Physiological noise modeling’, PNM (Brooks et al., 2008). Regressors were created for cardiac and respiratory signals using sine and cosine terms of their principal frequencies θC and θR, as well as their first three harmonics (2θC to 4θC and 2θR to 4θR). In addition, interactions of cardiac and respiratory signals were modeled using multiplicative terms of the form sin or cos (AθC ± BθC) with A,B = 1,2 (Kong et al., 2012). Altogether, 32 nuisance regressors were included and treated as voxel-wise confounds in FSL FEAT. As PNM is tailored for application in a general linear model, we had to choose a regressor of interest for the regression step. As low frequencies are of no importance in fMRI experiments, we chose a linear ramp function for this. Masked independent component analysis (mICA) Pre-processed fMRI data of the discovery sample (n = 50) were temporally concatenated and analyzed by probabilistic independent component analysis, pICA (Beckmann and Smith, 2004) using MELODIC 3.13 (FSL5.0). To suppress physiological noise in the spatial domain, pICA was restricted to the volume inside an anatomical brainstem mask, thus, excluding the most problematic areas adjacent to the brainstem from the analysis (masked ICA, mICA, Fig. 1f). The brainstem mask was defined based on gray and white matter tissue maps supplied with the ICBM152 template. The tissue probability was thresholded at 0.9 and remaining non-brainstem regions were removed manually (Supplementary Fig. 1). The relatively high threshold was chosen to result in a smaller mask in order to remove as many CSF-containing (and therefore noisy) voxels from the mask as possible. The final mask comprised the entire medulla, pons, and mesencephalon. Data were projected into a 37-dimensional subspace using probabilistic principal component analysis after voxel-wise de-meaning of the
93
data and normalization of the voxel-wise variance. The decomposition dimensionality was based on the median of all single subject dimensionality estimates using the Laplace approximation to the Bayesian evidence of the model order (Beckmann and Smith, 2004; Minka, 2000). Note that we did not use dimensionality estimates from the whole sample of n = 50, as it has been shown that estimates from analytic methods are strongly influenced by sample size (a higher number of observations leading to higher dimensionality estimates, although the underlying dimensionality stays the same) (Cordes and Nandy, 2006). The whitened observations were decomposed into time-courses and maps by optimizing for non-Gaussian spatial distributions using a fixedpoint iteration technique (Hyvärinen, 1999). Estimated component maps were divided by the standard deviation of the residual noise and thresholded by fitting a mixture model to the histogram of intensity (Beckmann and Smith, 2004). The analysis was repeated with the data of the confirmation sample (n = 50). Comparing mICA to existing methods To illustrate the performance of mICA compared to standard ICA in combination with other noise suppression methods, four separate ICAs were run on the cuboid containing the brainstem (Fig. 1): (1) Unmasked pICA of data without correction for physiological noise, (2) unmasked pICA after regression of physiological signals using PNM, (3) unmasked pICA after low-pass filtering with a 0.1 Hz cutoff, and (4) mICA after regression of physiological signals using PNM. Dimensionalities were estimated on the group level for each analysis using the Laplace approximation to the Bayesian evidence of the model order. The three independent components with the highest uniquely explained variance in each approach were extracted for visual comparison. Furthermore, variance maps were calculated from each of the four functional datasets. Reproducibility and symmetry analysis of mICA results Reproducibility analysis of the brainstem ICN was accomplished by calculating spatial cross-correlation coefficients for all possible combinations of the un-thresholded independent components from the discovery and the confirmation sample. The resulting cross-correlation matrix was then sorted using the Hungarian sorting algorithm for a linear assignment problem (Munkres, 1957). This yielded a one-to-one mapping maximizing the sum of all coefficients. Components with a correlation coefficient above 0.5 were considered reproducible. Symmetry of the brainstem ICN in the discovery sample was assessed following an approach similar to reproducibility analysis. Correlation coefficients were calculated for all components of the discovery sample and their left–right flipped counterparts. Again, Hungarian sorting yielded a one-to-one mapping, where components with an intrinsic symmetry were mapped on themselves and all other components on their best-matching symmetric partners. Components with a correlation coefficient larger than 0.5 were considered intrinsically symmetric and symmetric partners, respectively. Back-projection analysis to study cortical connectivity and specificity of brainstem components We back-projected the associated time-courses of all 37 mICA components of the discovery sample onto the whole-brain functional data (Calhoun et al., 2001) using a general linear model (Fig. 2). This was done to find all brain areas exhibiting a temporal behavior similar to that of the brainstem components under consideration. The reason for this was twofold. Firstly, we wanted to differentiate specific components (i.e. those of neuronal origin) from unspecific ones (i.e. those related to noise) (Fig. 2). Secondly, we were interested in the cortical and subcortical functional connectivity profile of the brainstem
94
F. Beissner et al. / NeuroImage 86 (2014) 91–98
proportions. Specificity was defined as s = 1-pCSF. Components were considered specific, if their s-value was within one standard deviation of the mean s. mICA vs. mICA + PNM To assess the impact of physiological noise regression (PNM) as a pre-processing step, when using the masked ICA approach, we also repeated mICA without prior noise regression for both discovery and confirmation sample. For the sake of comparability, the number of components was kept at 37. Matching components between the two samples (with and without PNM) were once more identified by calculating spatial cross-correlation coefficients and subsequent Hungarian sorting. We assessed reproducibility of components derived without prior PNM using the same approach as described for the PNM sample. A two-sided two-sample t-test was calculated for the hypothesis of reproducibility differences between the two analyses. A p-value of 0.05 was considered significant. Results mICA compared to existing methods
Fig. 2. Conceptual overview of the approach used to assess specificity of the brainstem components (see Fig. 4 for results). For each component, the associated time-course was back-projected into the whole-brain data to identify regions exhibiting a similar timecourse. The sum of these regions can be interpreted as the cortical (and sub-cortical) connectivity profile of the brainstem component. In a second step the thresholded connectivity patterns were scalarly multiplied with a probability mask of cerebro-spinal fluid (CSF) to calculate the percentage of activated voxels that lay within CSF regions (pCSF). Specificity was defined as s = 1-pCSF.
components. The analysis was done separately for each subject by splitting the concatenated time-courses of the ICA into single-subject timecourses. The single-subject parameter estimates were then tested for statistically significant areas using a non-parametric permutation test implemented in RANDOMISE 2.9 (FSL5.0). 5000 permutations of the data were calculated. A voxel-wise p-value of 0.01, corrected for multiple comparisons by controlling the family-wise error (FWE), was considered significant. For specificity analysis, we calculated the percentage pCSF of the overall volume of activated voxels that fell into cerebro-spinal fluid (CSF) regions as compared to white and gray matter. The rationale for this was that a high value would indicate a strong contribution of CSF signal to the brainstem component, supporting a non-neuronal origin. The thresholded, binarized p-maps were scalarly multiplied with the tissue probability map of CSF distribution available for the ICBM152 template. The same procedure was repeated for white and gray matter maps. Prior to multiplication all probability maps had been thresholded at a tissue probability of 0.5. pCSF was defined as the CSF proportion of activated voxels normalized by the sum of the white and gray matter
The results of the analysis under the three different conditions (without physiological noise correction, after regression of physiological signals, and after temporal low-pass filtering) are shown in Fig. 1. Dimensionality estimates were 16 for no correction, 22 for physiological noise regression, and 55 for low-pass filtering. In contrast, the dimensionality estimate for the mICA approach was 95. For each analysis, the three independent components with the highest uniquely explained variance are shown, while the entire set of components can be found in Supplementary Figs. 2, 3, and 4 (for the alternative approaches). The dominant independent components of all three alternative analyses showed similar artifacts, which could easily be sub-divided into two classes. The majority of ICA components showed a pattern of activity that was chiefly located outside the brainstem. A comparison with physiological noise patterns (Fig. 1b) revealed a strong similarity between their spatial distribution and that of the ICA results. This is expected, as ICA algorithms identify signal sources by their nonGaussianity, a property shared by almost all physiological noise signals. As noise signals have stronger signal intensities than BOLD signals, they dominate the ICA results. The second class of artifacts was also found in the results of all three analyses: Extensive activation clusters that were partly located in the brainstem, but primarily comprised subcortical and cerebellar regions. These components result from BOLD signals in nonbrainstem regions whose signal intensities exceed those within the brainstem. Although they are valid components depicting neuronal activity, they cannot be used to study the behavior of brainstem nuclei, as they usually form large clusters (cf. Supplementary Fig. 3, lower left corner). When applying an anatomical brainstem mask, and restricting the ICA within that mask (i.e. the mICA approach), the analysis yielded starkly contrasting results: as shown in Fig. 1f, activation clusters of single nuclei and nuclear complexes were made evident, and a relative lack of artifacts could be observed (Fig. 1f). Reproducibility and symmetry of mICA results The results of the mICA analysis for the discovery sample are shown in Fig. 3. For comparison, those of the confirmation sample can be found in Supplementary Fig. 5. Components were ordered by their anatomical localization from rostral to caudal and show a rather equal distribution between mesencephalon (10 ICs), pons (20 ICs), and medulla (7 ICs) taking into account the relative volume of these parts of the brainstem.
F. Beissner et al. / NeuroImage 86 (2014) 91–98
95
Fig. 3. Results of a 37-dimensional masked independent component analysis of 50 healthy subjects. Brainstem components are presented at a mixture model threshold of 0.5 and ordered in rostro-caudal direction. Reproducibility (rep) was assessed by spatial cross-correlation with components in a confirmation sample (Supplementary Fig. 5), while specificity (spec) scores are based on whole-brain connectivity of the components (Fig. 2). Graphs on the lower right show the distributions of these scores as well as the cutoff used to identify non-reproducible (ο) and non-specific (i.e. noise-related) components (#). Note that the ordering of the components is different in both graphs.
The majority (32/37) of all components in the discovery sample were reproducible, i.e. they had a correlation coefficient higher than r = 0.5 with their best matching partner from the confirmation sample. The average correlation coefficient of the reproducible components was 0.65 ± 0.11 (mean ± STD). Components showed a high symmetry, in that 32 of the components were either intrinsically symmetric (14 components) or had a contralateral partner (18 components). This means that all these components had a correlation coefficient higher than r = 0.5 with one of the other components of the sample. The average correlation coefficient of the symmetric components was 0.64 ± 0.08 (mean ± STD).
Whole-brain connectivity and specificity of brainstem components Whole-brain connectivity profiles showed a large variability across brainstem components. A group of representative components is shown in Fig. 4. Some of them resembled known large-scale brain networks, like the default mode network (Fig. 4a), while others did not (Figs. 4b+c). Three classes of components could be distinguished by eye: The first contained the abovementioned networks of cortical and subcortical structures, where connectivity was clearly dominated by gray matter regions (Figs. 4a–c). In the second class, profiles comprised a much higher percentage of CSF and white matter regions (Figs. 4d–f).
96
F. Beissner et al. / NeuroImage 86 (2014) 91–98
Fig. 4. Results of the cortical connectivity/specificity analysis (see Fig. 2 for method). Three representative components are shown for the specific and unspecific group together with their specificity values (spec). While the cortical connectivity profiles of the specific components (a–c) are dominated by gray matter regions, unspecific components (d–f) exhibit noticeable connectivity with CSF regions. An interesting observation that further corroborates the validity of our approach is the median pontine component in (a) showing a connectivity profile that largely resembles the default mode network, a connection that has been reported before (Habas et al., 2009).
Specificity analysis confirmed this distinction, as the latter components had higher pCSF and therewith lower specificity values. The third class comprised only two components, namely the caudal and rostral periaqueductal gray, which showed only very localized connectivity that was, however, dominated by CSF of the aqueduct leading to their classification as non-specific. mICA vs. mICA + PNM Omitting physiological noise regression in the mICA analysis yielded results that were very similar to that of the original analysis (Supplementary Fig. 6). The average correlation coefficient between both sets of components was 0.79 ± 0.16 (mean ± STD). Only two components had correlation coefficients below 0.5. One of them was a bilateral activation pattern of the superior colliculi, the other one a putative noise component. Reproducibility did not significantly differ between the PNM and non-PNM analyses (0.59 ± 0.21 vs. 0.61 ± 0.21, p = 0.46) (Supplementary Fig. 6). Discussion In this study, we combined a multivariate data-driven analysis method (ICA) with an anatomical mask and applied it to resting-state fMRI data, to study the activity of brainstem nuclei and intrinsic connectivity networks. To the best of our knowledge, it constitutes the first successful attempt to depict brainstem nuclei as well as their connectivity to cortical areas in humans using ICA. The key distinction between our masked ICA approach and previous methods is the radically different treatment of physiological noise by exploiting mainly spatial instead of temporal characteristics.
Instead of using a two-step procedure of mICA and back-projection to arrive at brainstem–cortex connectivity profiles, one could also think of using a whole-brain ICA approach. There are, however, several reasons, why this common approach usually does not identify brainstem foci associated with the cortical component. One is that whole-brain group analyses require at least a minimal amount of spatial smoothing. This pre-processing step inevitably leads to contamination of the brainstem with physiological noise from the adjacent CSFcontaining regions. Furthermore, the optimal amount of smoothing to depict brainstem nuclei may be different from that commonly used for cortical studies (Beissner et al., 2011). The most important reason, why brainstem- and whole-brain-centered approaches produce different results, is that whole-brain ICAs are largely driven by cortical signals due to the much larger relative volume of the cortex compared to the brainstem. Thus, brainstem nuclei are only detected, if they belong to one of the known large-scale resting state networks. Fig. 4a shows such an example, where back-projection of a median pontine component yielded the default mode network — a result that corroborates previous reports of pontine nuclei being part of the default mode network (Habas et al., 2009). In general, however, the cortical counterpart of brainstem networks may be small compared to large-scale brain networks and may, thus, remain completely undetected by whole-brain ICA. As Figs. 4b+c show, such cortical components may indeed look very different from those commonly reported in the literature. Concerning reproducibility and specificity of our results, we found that brainstem resting-state components show high reproducibility, at least across samples acquired at the same MR scanner. However, our results also show that reproducibility must not be equated with specificity. The fact that the five components with the lowest specificity scores all showed moderate reproducibility, underlines that even components
F. Beissner et al. / NeuroImage 86 (2014) 91–98
that are reproducible may still be of non-neuronal, i.e. artifactual, origin. This result, however, is not unexpected, as physiological noise is highly structured and there is no reason to believe that its spatial profile should differ between subjects. As mICA relies mainly on a spatial approach to remove physiological noise (by excluding highly affected areas from the analysis), we have combined it with PNM, a temporal noise reduction method (Brooks et al., 2008). As the two approaches are complementary, one could think that their combination should yield better results than either of it alone. However, while pre-processing with PNM led to a higher dimensionality estimate compared to the uncorrected analysis (Supplementary Figs. 2 vs. 3), we were not able to find significant differences in reproducibility when a fixed dimensionality was used for the comparison (Supplementary Fig. 6). The decision whether or not an additional temporal noise regression step should be applied to the data has to be made individually for each study. On the one hand, a combination of both techniques (PNM + mICA) may be a more conservative approach, whose potential should be assessed in further studies. On the other hand, every regression of physiological “noise” bears the risk of removing meaningful signal from the data (Iacovella and Hasson, 2011; Khalili-Mahani et al., 2012). This is especially relevant for studies of autonomic nuclei or of those triggering autonomic responses (e.g. nociceptive). One great advantage of omitting PNM is of course the applicability of mICA to existing fMRI data sets that lack physiological recordings. However, the authors should take great care to control for possible residual noise. Another pre-processing step that may be optional is the upsampling of the brainstem data to 1 mm. While this step clearly improves post-hoc comparison of functional and structural data, a resolution of 2 mm was found to produce very similar results later (results not shown). Several limitations should be addressed. Firstly, we used a restingstate experiment, whereas specific tasks may have facilitated the identification of nuclei. However, correspondence of the brain's functional architecture during tasks and resting-state has recently been demonstrated for the cortex (Smith et al., 2009), a concept that can probably be transferred to the brainstem. Furthermore, validation of brainstem nuclei based on task-fMRI would be hampered by the lack of tasks that activate only a small number of nuclei. For example, letting subjects move their eyes to a fixed cue would not only activate the abducens nuclei, but also visual, attention, vestibular and many more systems localized side-by-side in the brainstem. Another problem is that almost all nuclei have multiple functions (Koutcherov et al., 2009) and for many of them, these functions are still unknown. Secondly, by restricting the ICA to the brainstem region, we lost spatial information from the rest of the brain that could otherwise be used to identify imaging artifacts (e.g., ghosting and motion). However, parts of the information lost by cropping of the images are regained by our back-projection approach. As Supplementary Figs. 2–4 show, artifacts in the brainstem region also have specific patterns making their visual identification possible. Our results raise the prospect that the measurement of intrinsic and extrinsic connectivity of the brainstem will become standard in fMRI studies, both in humans and in small animal imaging and complement information from cortical and subcortical regions. Studies of brainstem–cortex connectivity in humans will prove crucial for the understanding of functional systems, like the neuromodulatory (Briand et al., 2007) and nociceptive (Heinricher et al., 2009) systems or the central autonomic network (Benarroch, 1993), as all of them have vital nodes in the brainstem, and extensive connections with the cortex. Furthermore, these systems have critical clinical relevance; notable examples include the neuromodulatory systems in major depression (Belmaker and Agam, 2008) and Parkinson's disease (Obeso et al., 2010), the nociceptive system in chronic pain syndromes (Tracey and Mantyh, 2007), and the autonomic nervous system in hypertension and cardiovascular diseases (Brook and Julius, 2000).
97
Finally, it is the hope of the authors that the brainstem will lose its label of a terra incognita and soon become a region of major interest in the neuroimaging community. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2013.07.081.
Acknowledgments FBe thanks Vitaly Napadow for helpful discussions. FBe was supported by the German Research Foundation Grant BE4677/1-1. FBr and AS were supported by the Horst-Görtz-Foundation.
References Ashburner, J., Friston, K.J., 2005. Unified segmentation. NeuroImage 26 (3), 839–851. http://dx.doi.org/10.1016/j.neuroimage.2005.02.018. Beckmann, C.F., Smith, S.M., 2004. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23 (2), 137–152. Beissner, F., Deichmann, R., Baudrexel, S., 2011. fMRI of the brainstem using dual-echo EPI. NeuroImage 55 (4), 1593–1599. http://dx.doi.org/10.1016/j.neuroimage.2011.01.042. Belmaker, R.H., Agam, G., 2008. Major depressive disorder. N. Engl. J. Med. 358 (1), 55–68. http://dx.doi.org/10.1056/NEJMra073096. Benarroch, E., 1993. The central autonomic network: functional organization, dysfunction, and perspective. Mayo Clin. Proc. Mayo Clin. 68, 988. Briand, L.A., Gritton, H., Howe, W.M., Young, D.A., Sarter, M., 2007. Modulators in concert for cognition: modulator interactions in the prefrontal cortex. Prog. Neurobiol. 83 (2), 69–91. http://dx.doi.org/10.1016/j.pneurobio.2007.06.007. Brook, R.D., Julius, S., 2000. Autonomic imbalance, hypertension, and cardiovascular risk. Am. J. Hypertens. 13 (6 Pt 2), 112S–122S. Brooks, J.C.W., Beckmann, C.F., Miller, K.L., Wise, R.G., Porro, C.A., Tracey, I., Jenkinson, M., 2008. Physiological noise modelling for spinal functional magnetic resonance imaging studies. NeuroImage 39 (2), 680–692. http://dx.doi.org/10.1016/j.neuroimage. 2007.09.018. Calhoun, V.D., Adali, T., Pearlson, G.D., Pekar, J.J., 2001. A method for making group inferences from functional MRI data using independent component analysis. Hum. Brain 151, 140–151. http://dx.doi.org/10.1002/hbm. Cordes, D., Nandy, R., 2006. Estimation of the intrinsic dimensionality of fMRI data. NeuroImage 29 (1), 145–154. D'Ardenne, K., McClure, S.M., Nystrom, L.E., Cohen, J.D., 2008. BOLD responses reflecting dopaminergic signals in the human ventral tegmental area. Science 319 (5867), 1264–1267. http://dx.doi.org/10.1126/science.1150605. Eippert, F., Bingel, U., Schoell, E.D., Yacubian, J., Klinger, R., Lorenz, J., Büchel, C., 2009. Activation of the opioidergic descending pain control system underlies placebo analgesia. Neuron 63 (4), 533–543. http://dx.doi.org/10.1016/j.neuron.2009.07.014. Formisano, E., Esposito, F., Di Salle, F., Goebel, R., 2004. Cortex-based independent component analysis of fMRI time series. Magn. Reson. Imaging 22 (10), 1493–1504. http:// dx.doi.org/10.1016/j.mri.2004.10.020. Greve, D.N., Fischl, B., 2009. Accurate and robust brain image alignment using boundarybased registration. NeuroImage 48 (1), 63–72. http://dx.doi.org/10.1016/j.neuroimage. 2009.06.060. Griswold, M.A., Jakob, P.M., Heidemann, R.M., Nittka, M., Jellus, V., Wang, J., Kiefer, B., et al., 2002. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn. Reson. Med. 47 (6), 1202–1210. http://dx.doi.org/10.1002/mrm.10171. Habas, C., Kamdar, N., Nguyen, D., Prater, K., Beckmann, C.F., Menon, V., Greicius, M.D., 2009. Distinct cerebellar contributions to intrinsic connectivity networks. J. Neurosci. 29 (26), 8586–8594. http://dx.doi.org/10.1523/JNEUROSCI.1868-09.2009. Harvey, A.K., Pattinson, K.T.S., Brooks, J.C.W., Mayhew, S.D., Jenkinson, M., Wise, R.G., 2008. Brainstem functional magnetic resonance imaging: disentangling signal from physiological noise. J. Magn. Reson. Imaging 28 (6), 1337–1344. Heinricher, M.M., Tavares, I., Leith, J.L., Lumb, B.M., 2009. Descending control of nociception: specificity, recruitment and plasticity. Brain Res. Rev. 60 (1), 214–225. http://dx.doi.org/10.1016/j.brainresrev.2008.12.009. Hyvärinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Netw. 10 (3), 626–634. Iacovella, V., Hasson, U., 2011. The relationship between BOLD signal and autonomic nervous system functions: implications for processing of “physiological noise”. Magn. Reson. Imaging 29, 1338–1345. Khalili-Mahani, N., Chang, C., Van Osch, M.J., Veer, I.M., Van Buchem, M.A., Dahan, A., Beckmann, C.F., et al., 2012. The impact of “physiological correction” on functional connectivity analysis of pharmacological resting state fMRI. NeuroImage 65, 499–510. http://dx.doi.org/10.1016/j.neuroimage.2012.09.044. Klose, U., Strik, C., Kiefer, C., Grodd, W., 2000. Detection of a relation between respiration and CSF pulsation with an echoplanar technique. J. Magn. Reson. Imaging:JMRI2 11 (4), 438–444. Kong, Y., Jenkinson, M., Andersson, J., Tracey, I., Brooks, J.C.W., 2012. Assessment of physiological noise modelling methods for functional imaging of the spinal cord. NeuroImage 60 (2), 1538–1549. http://dx.doi.org/10.1016/j.neuroimage.2011.11.077. Koutcherov, Y., Huang, X.F., Halliday, G., Paxinos, G.T., 2009. Organization of human brain stem nuclei. In: Paxinos, G.T., Mai, J.K. (Eds.), The Human Nervous System. Elsevier, London, pp. 267–320.
98
F. Beissner et al. / NeuroImage 86 (2014) 91–98
McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S., Bell, A.J., Sejnowski, T.J., 1998. Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapp. 6 (3), 160–188. Minka, T., 2000. Automatic choice of dimensionality for PCA. Technical Report 514. Munkres, J., 1957. Algorithms for the assignment and transportation problems. J. Soc. Ind. Appl. Math. 5 (1), 32–38. Obeso, J.A., Rodriguez-Oroz, M.C., Goetz, C.G., Marin, C., Kordower, J.H., Rodriguez, M., Hirsch, E.C., et al., 2010. Missing pieces in the Parkinson's disease puzzle. Nat. Med. 16 (6), 653–661. http://dx.doi.org/10.1038/nm.2165. Smith, S.M., Fox, P., Miller, K., 2009. Correspondence of the brain's functional architecture during activation and rest. Proceedings of the National Academy of Sciences of
the United States of America 106 (31), 13040–13045. http://dx.doi.org/10.1073/ pnas.0905267106. Thompson, S.K., Von Kriegstein, K., Deane-Pratt, A., Marquardt, T., Deichmann, R., Griffiths, T.D., McAlpine, D., 2006. Representation of interaural time delay in the human auditory midbrain. Nat. Neurosci. 9 (9), 1096–1098. http://dx.doi.org/10.1038/ nn1755. Tracey, I., Mantyh, P.W., 2007. The cerebral signature for pain perception and its modulation. Neuron 55 (3), 377–391. http://dx.doi.org/10.1016/j.neuron.2007.07. 012. Wessel, N., 2000. Nonlinear analysis of complex phenomena in cardiological data. Herzschrittmacherther. Elektrophysiol. 11, 159–173.