The Subcortical Cocktail Problem; Mixed Signals from the Subthalamic Nucleus and Substantia Nigra

The subthalamic nucleus and the directly adjacent substantia nigra are small and important structures in the basal ganglia. Functional magnetic resonance imaging studies have shown that the subthalamic nucleus and substantia nigra are selectively involved in response inhibition, conflict processing, and adjusting global and selective response thresholds. However, imaging these nuclei is complex, because they are in such close proximity, they can vary in location, and are very small relative to the resolution of most fMRI sequences. Here, we investigated the consistency in localization of these nuclei in BOLD fMRI studies, comparing reported coordinates with probabilistic atlas maps of young human participants derived from ultra-high resolution 7T MRI scanning. We show that the fMRI signal reported in previous studies is likely not unequivocally arising from the subthalamic nucleus but represents a mixture of subthalamic nucleus, substantia nigra, and surrounding tissue. Using a simulation study, we also tested to what extent spatial smoothing, often used in fMRI preprocessing pipelines, influences the mixture of BOLD signals. We propose concrete steps how to analyze fMRI BOLD data to allow inferences about the functional role of small subcortical nuclei like the subthalamic nucleus and substantia nigra.


Introduction
The subthalamic nucleus (STN) and substantia nigra (SN) are small subcortical nuclei in the basal ganglia [1]. Many cognitive neuroscience studies using BOLD functional magnetic resonance imaging (fMRI) have shown that the STN and SN are involved in a range of tasks such as response inhibition [2][3][4], conflict processing [5], force production [6,7], working memory [8,9], and the adjustment of global and selective response thresholds [10] (Table 1). While several studies have shown the involvement of either the STN or the SN in tasks such as the stopsignal task or Simon task, several neurocomputational models make distinct functional predictions for these nuclei (e.g. [11,12]). For instance, Frank et al. [12,13] propose that the STN acts as a general brake, whereas the SN, through the release of dopamine, is activated by the correct response and inhibits the incorrect response. However, the small size and close proximity of the STN to the SN, as well as to the surrounding brain structures, frustrates precise localization and makes it challenging to attribute the BOLD signal to either STN or SN, respectively. A common procedure is to place a box of 10x10x10 mm on a center coordinate in the STN and extract the mean BOLD fMRI signal in this box, to estimate the signal change in the STN [2,14,15]. Note that the volume of such a box is about 7.5 times larger than the average STN volume reported in the literature (1000 mm 3 as compared to a weighted average of 119.88 mm 3 and a weighted median of 131.75 mm 3 [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. In combination with often-used smoothing procedures, the signals originating from the STN and SN will get mixed, making it difficult to unequivocally attribute signal to either structure [30,31]. The goal of this study was to investigate the consistency of the coordinates found in fMRI studies on the STN and SN, summarize the methods employed in these studies, and assess the severity of the problems with localization and mixture of signals. In a first step, we conducted a comprehensive literature search to characterize the methods resulting in significant functional activation in the STN and SN. In a second step, the peak coordinates of the STN and SN derived from these studies were compared to the location of recently published probability STN and SN ultra-high resolution 7T MRI atlas [28]. Thirdly, using ultra-high resolution individual anatomical MRI masks, a simulation study was performed to test the influence of different smoothing kernels on the mixture of BOLD fMRI signals from both the STN and SN.

Selection of STN and SN BOLD fMRI studies
A comprehensive search for relevant neuroimaging studies in the field of BOLD fMRI studies including the STN and SN was carried out using Google scholar (http://scholar.google.com/). The main keywords utilized were 'fMRI + substantia nigra', 'fMRI + SN', 'fMRI + subthalamic nucleus', 'fMRI + STN', as well as all combinations of the aforementioned terms.
Based on the information contained in the abstracts of all the papers returned, empirical studies were selected to meet the following inclusion criteria: (1) Studies were published in peer-review English language journals between January 2000 and March 2014; (2) the studies used BOLD fMRI; (3) the studies reported a functional coordinate that could be attributed to either the SN or STN; and (4) the studies reported the location of activation as 3D coordinates in stereotactic space of Talairach or the Montreal Neurological Institute (MNI).
All empirical studies included were cross-referenced and all papers citing these empirical studies were searched, using the Google scholar citation index tool. The whole selection process was repeated for the newly obtained empirical papers until no new studies were found. This resulted in the inclusion of 52 papers (see Table 1).
All activation foci of the included studies that were originally reported in Talairach space were converted to the MNI stereotactic space using the Lancaster et al. transformation algorithm, which has been validated and shown to substantially reduce any bias between the two references spaces [32].

Probabilistic ultra-high resolution 7T MRI atlas maps
For analysis of the comparison between STN and SN coordinates reported in the literature (see Table 1), we used previously reported ultra-high resolution 7T MRI probability maps [28]. The probability maps are based on 30 participants (14 females) with a mean age of 24.2 year (SD 2.4). The STN and SN masks were manually segmented by two raters for each individual on 7T zoomed multi-echo 3D FLASH MRI data with an isotropic voxel size of 0.5 mm [33]. Only voxels rated by both raters as belonging to the STN or SN were included in further analyses. Note that no differentiation between the SN pars compacta and the pars reticulate were made because the voxel resolution and used scan sequence did not allow for identification of the two subparts. The individual masks were then linearly registered to MNI standard space and combined to create a probabilistic atlas. For more information regarding the segmentation, MRI scanning sequence, and registration procedure see [28,34]. The structural data can be found on http://www.nitrc.org/projects/atag_mri_scans/ and on http://dx.doi.org/10.5061/dryad.fb41s. The probabilistic masks can be found on http://www.nitrc.org/projects/atag.

Simulation study
A simulation study was performed to assess the amount of signal that originates from neighboring nuclei that can be introduced into a region of interest (ROI) by smoothing. Sixty STN and sixty SN masks (thirty masks in both hemispheres) from the ATAG (Atlas of The bAsal Ganglia) dataset [28] were used in a total of 60 simulations, all using one STN and one SN mask at a voxel resolution of 0.5 mm isotropic. It was assumed that every voxel in each mask contained a signal of unit strength. Then, smoothing kernels of different sizes were applied, and for every voxel and for every nucleus, the amount of signal in that voxel originating from that nucleus was determined. If the entire signal came from the same nucleus, the value was 1. If no signal from that nucleus reached that voxel, the value was 0. The sum of the signal strengths of the two nuclei in a voxel could never surpass 1.
We focused on the mixture of signal in the center voxel of both masks, to emulate a ROI study where the ROI would be placed in the best possible voxel according to the ground truth. This is a very optimistic scenario considering the difficulty of STN/SN localization as discussed earlier. For every center voxel, two quantities were computed: corresponding to the amount of signal originating from the SN and the amount of signal originating from the STN.
[x com ,y com ,z com ] is the coordinate of the center-of-mass of the mask-of-interest in millimeters.STN_mask[x,y,z] was either 1 or 0, corresponding to the coordinate [x,y,z] being in the STN or not, SN_mask[x,y,z] analogously for the SN. σ is the standard deviation of the Gaussian kernel, which can be calculated for a given FWHM (full width at half maximum) by using the following formula:

Overview of functional MRI STN and SN studies
52 functional MRI studies were included in the present study (Table 1), published between 2003 and 2014. These studies employed for instance resting state, the stop-signal task, decision-making tasks including reward and outcome manipulations, and threshold adjustments in cognitive control tasks. Smoothing kernels ranged between 0 mm -12 mm and data was collected on either 1.5 or 3T scanners. Note that not all studies report extensive methodological or procedural information, which limits the assessment of their anatomical specificity [35].

STN and SN coordinates in MNI space compared to ultra-high resolution 7T atlas maps
The reported STN and SN-coordinates were compared to the center-of-mass coordinates of the previously published 7T MRI probabilistic masks [28]. Results are summarized in Table 2.
On average, reported STN activity coordinates lay 5.2 mm (left hemisphere, std. = 3.3) and 5.7 mm (right hemisphere, std. = 2.8) from the center-of-mass of the 7T MRI probabilistic mask. Several studies include older participants [15,[36][37][38]). This might result in a mismatch between the reported coordinates and the probabilistic atlas because it is known that the STN shifts in lateral direction with age [29,39,40]. Note, however, that this lateral shift is smaller (on average 1.6 mm more lateral for elderly than for young participants [29]) compared to the standard deviation reported in the present study. The mismatch in location is predominantly observed in the dorsal-ventral and anterior-posterior direction such that the reported coordinates were on average 1.3 / 1.5 (left/right, std. = 2.8/2.5) mm more dorsal and 1.2 / 2.5 (left/right std. = 4.4/4.0) mm more anterior than the center-of-mass of the 7T probabilistic STN masks.
The left STN coordinates were on average 0.9 mm (std. = 2.8) more medial than the centerof-mask of the 7T MRI probabilistic mask, and the right STN coordinates lay on average 0.4 mm (std. = 3.2) more lateral. The average distance between the centers-of-mass of the SN and STN of the probabilistic masks is in the same order of magnitude as the distance between the average reported STN fMRI location and the actual STN center-of-mass (6.4 mm left, 6.7 mm right; std. = 0.7 / 0.7).
Reported SN activity coordinates were on average 5.6 (left hemisphere, std. = 4.0) and 8.9 (right hemisphere, std. = 8.6) mm away from the center of mass of the 7T MRI probabilistic mask. Shifts occurred in all three directions. Fig. 1 shows coronal plots of the probabilistic maps with the reported coordinates rendered onto them. Only coordinates that lay within the MNI coordinate range of x: 18,-18 / z: 0,-21 / y: −2,-25 were plotted. The coordinates that fell outside of this range are marked in Table 1.
The average FWHM size of the smoothing kernel was 6.3 mm (median 6 mm). The 8 mm FWHM smoothing kernel was used most frequently (16 out of 52 studies: see Fig. 2 for the relative size of the kernels used compared to the STN and SN). There was no relationship between the nucleus of interest and the smoothing kernel used. To test for any spatial biases introduced by smoothing, as reported by e.g. Sacchet et al., [41] for the nucleus accumbens, MNI coordinates were correlated with the size of the smoothing kernel employed. No correlations were found except for the SN (r(45) = .83, p <0.05, uncorrected), which lay more superior as a function of a larger smoothing kernel. The majority of studies (47 out of the 52) reported the voxel resolution. The voxel resolution was on average 34 mm 3 (median 32 mm 3 , std. = 21 mm, range 3.4-103.8 mm). The result shows that signal originating from each nucleus spreads widely, also across its neighbor's boundaries (Fig. 4).  Table 1. Overlaid onto these coordinates is the probabilistic atlas of the STN and SN. The isolines reflect the percentage overlap across the 30 young subjects taken from Keuken et al. [28]. The outermost isolines reflects a 10% probability of containing the SN at the population level, the more inner lines represent 30%, 50%, and 70% probability of containing the SN. The outermost isolines for the STN reflects a probability of 20% containing the STN, the inner line represent 40% probability of containing the STN. The grid size corresponds to a voxel size of 3x3 mm. All coordinates are in MNI standard space. When no smoothing was applied, the signal in the center voxel of both SN and STN originated completely from the nucleus of its location. When a 4 mm FWHM smoothing kernel was used, 30% of the signal in the center voxel of the STN was found to originate from outside the STN and SN, while ten percent originated from the SN and 60% from the STN itself. With an often-used 8 mm FWHM smoothing kernel, 75% of the signal in the center voxel of the STN mask originated from outside the STN and SN. Ten percent originated from the SN and only 15% from the STN itself. For this simulation, the strength of the STN signal in the center voxel was taken to be similar in size to the signal originating from the SN. Note that the results would vary if different signal strengths in both nuclei were assumed. However, the precise ratios have little effect on this finding. Importantly, this simulation shows that with large smoothing kernels it becomes impossible to disentangle the origin of the measured signal, even when focusing on the most central voxel. In empirical fMRI studies this is an unlikely scenario, because the voxel resolution is at an average of 34 mm3, instead of the 0.064 mm3 (0.4 mm isotropic) used here. In sum, the mixing of signals is likely to be substantially worse in empirical fMRI studies compared to this simulation study.

Discussion
In the present study we show that there is large variability in previously reported fMRI coordinates attributed to the STN and SN. We also show a discrepancy between individual coordinates of empirical studies and probabilistic atlas maps derived from ultra-high resolution 7T MRI [28]. The resolution of the fMRI sequences used in the studies we included was usually  The simulation results reveal that when smoothing kernels of commonly used sizes are applied, the amount of signal from neighboring nuclei that get smoothed into a region of interest is of similar size as the signal from the region itself. This is particularly important when analyzing data from a small nucleus such as the STN, which borders the larger SN.
These results add empirical data to the recent discussion about smoothing in functional neuroimaging. Stelzer et al. [31] suggested that smoothing fMRI data should be abandoned altogether, because it (1) causes incorrect estimation of the true spatial extent of brain activation, (2) blurs away signals of limited spatial extent, and (3) frustrates the detection of low-intensity signals in the vicinity of non-active tissue. Our results illustrate quantitatively how large these effects can be, specifically for subcortical nuclei: we show that reported MNI coordinates largely non-overlap with anatomical masks (point 1) and that smoothing can induce substantial mixing with signal from outside the nucleus (point 2 and 3).
The use of smoothing can increase the signal-to-noise ratio in fMRI when the signal is more spatially correlated than the noise on the scale of the smoothing kernel employed. However, in the case of subcortical nuclei, the used smoothing kernels are often too large and mix in signal and noise from neighbouring structures. Yoon et al. [8] provide an empirical example of the influence of kernel size in their supplemental information: the activity in the SN only reached a significance threshold when a smoothing kernel with a very minor FWHM of 2 mm was employed. When a smoothing kernel with an FWHM of 8 mm was used, the effect disappeared. Because the voxel size of this study was rather large, 3.4x4x4 mm, the effect of a 2 mm  [28], the effect of smoothing on the mixing of signals in STN and SN in their respective center voxels was estimated. The lines show the amount of signal for different source-destination pairs of STN and SN as a function of smoothing kernel size. When no smoothing is applied, all signal in the SN originates from SN and all signal in STN originates from the STN. When more smoothing is applied, the amount of signal originating from the nucleus that is measured sharply decreases, and within the STN the amount of signal from the SN becomes equal in size to the signal originating from the STN itself.
doi:10.1371/journal.pone.0120572.g004 smoothing kernel was negligible and could have been abandoned altogether. When one applies such a relatively large smoothing kernel to data of such a relatively coarse resolution, the amount of signal in a voxel in the smoothed image originating from outside this voxel is less than 0.2% (See http://nbviewer.ipython.org/gist/Gilles86/0c093962de8cf05f76c8). The results by Yoon et al. thus clearly show that smoothing is not necessary to find significant effects in the substantia nigra region, even with 1.5T [8].
It has been suggested that a lack of spatial resolution and anatomical specificity could be overcome by using unsupervised clustering algorithms such as principal component analysis (PCA; [42]) or independent component analysis (ICA; [42]). These methods might 'detect' the nucleus of interest by exploiting the different covariance structures of the BOLD signal in different nuclei. We think, however, that such an approach is not appropriate. First and foremost, it assumes that the task-related BOLD activity in the STN and SN are uncorrelated. This is highly unlikely because both nuclei are part of the same functional networks, e.g., the basal ganglia motor control loops. Secondly, even if the signal could be separated to some extent, there is no objective way of finding out which cluster component belongs to which nucleus and to which extent they account for only one nucleus. Third, independent components might represent non-BOLD signals such as physiological noise. Fourth, the most adequate procedure of defining the actual signal of the nucleus of interest by means of, e.g. a demixing matrix (e.g., [43] or a Gaussian sphere [44]) remains elusive.
Therefore we suggest that during functional imaging of small subcortical nuclei, standard smoothing strategies should be avoided altogether. More complex, adaptive smoothing approaches [45] might be useful, but analysis protocols that do not require smoothing should be preferred. A-priori ROI analyses [46] do not require smoothing, nor do whole-brain univariate analysis approaches that make use of False Discovery Rate (FDR) as multiple comparison correction, as well as multivariate analysis strategies [31,47].
Concretely, we propose an approach that maximizes both anatomical specificity and signalto-noise. Researchers are advised to use individual anatomical masks based on an appropriate MR contrast (i.e., T2 Ã or quantitative susceptibility mapping (QSM)) that allows for detailed visibility and segmentation of the structures of interest (see, e.g., [27][28][29][48][49][50]). When individual segmentation is not feasible, researchers can use probabilistic atlas maps, as provided for the STN and SN in [27][28][29]. If the research question does not focus on anatomical patterns within the nucleus itself, the mean signal across all voxels in the nucleus can be analysed. This maximizes SNR and removes both the multiple comparisons problem, as well as the need for registration to a standard space. When different activation patterns within the nucleus are expected, a voxel-wise analysis within the anatomical mask can be computed.
Given the variability in reported coordinates and smoothing, one may question the validity of earlier fMRI findings in the STN/SN. It is important to note that studies with Parkinson Disease patients using deep-brain stimulation (DBS) [51] or lesioning of the STN [52] deliver important causal evidence for the functional role of the STN in motor control. fMRI studies that report BOLD activity in motor control paradigms are thus likely to be sensitive to actual task involvement. However, we believe that caution is warranted in interpreting the anatomical specificity of these findings. Especially interpreting findings from studies that 1) use smoothing kernels with a FWHM of more than 4 mm, 2) do not use anatomical masks that are based on individual anatomy, either individually segmented or based on a population probabilistic map like the ATAG dataset [28], and finally, 3) use voxel resolutions that are in the same order of magnitude as the nucleus itself.
In sum, the present study provides evidence for the importance of accounting for individual anatomy when attempting to understand the functional role of small subcortical areas such as the STN and SN. Moreover, the combination of ultra-high resolution fMRI with a very high voxel resolution and zoomed-in acquisition protocols will help to unmix signals arising from small subcortical structures in very close proximity. Finally, the simulation results indicate that spatial smoothing should be avoided when one values anatomical specificity of functional neuroimaging results.