Human Neuroscience

of the habenula in error monitoring and as a critical modulatory relay between the limbic forebrain structures and the midbrain. In another study, co-activation of the habenula and midbrain was observed during negative but not positive feedback (Shepard et al., 2006). A recent study with more detailed anatomical mapping confirmed the habenula as the locus responding to negative rewards (Salas et al., 2010). However, it remains unclear whether the habenula signals to the VTA/SN in outcome monitoring, as has been demonstrated in non-human primates. The current study aimed to fill this gap of knowledge. In previous studies, we observed activation of subcortical structures including the habenula during error trials in a stop signal task (Li et al., 2008b). Here we substantiated the functional con-nectivity of the habenula and VTA/SN, using psychophysiological interaction (PPI), Granger causality analysis (GCA), and mediation analysis. PPI is a voxel-wise method widely used to examine whether correlation in activity between two brain areas is modulated by psychological contexts (Friston et al. We used PPI to establish greater connectivity between the habenula and VTA/SN during stop error (SE) as compared to stop success (SS) trials. However, PPI does not specify the direction of influence between brain regions. In contrast, GCA has been used to model directional interaction between blood oxygenation level dependent

of the habenula in error monitoring and as a critical modulatory relay between the limbic forebrain structures and the midbrain. In another study, co-activation of the habenula and midbrain was observed during negative but not positive feedback (Shepard et al., 2006). A recent study with more detailed anatomical mapping confirmed the habenula as the locus responding to negative rewards (Salas et al., 2010). However, it remains unclear whether the habenula signals to the VTA/SN in outcome monitoring, as has been demonstrated in non-human primates. The current study aimed to fill this gap of knowledge.
In previous studies, we observed activation of subcortical structures including the habenula during error trials in a stop signal task (Li et al., 2008b). Here we substantiated the functional connectivity of the habenula and VTA/SN, using psychophysiological interaction (PPI), Granger causality analysis (GCA), and mediation analysis. PPI is a voxel-wise method widely used to examine whether correlation in activity between two brain areas is modulated by psychological contexts (Friston et al., 1997;Gitelman et al., 2003;Stephan et al., 2003;O'Reilly et al., 2008;Hare et al., 2009). We used PPI to establish greater connectivity between the habenula and VTA/SN during stop error (SE) as compared to stop success (SS) trials. However, PPI does not specify the direction of influence between brain regions. In contrast, GCA has been used to model directional interaction between blood oxygenation level dependent (BOLD) time series (Goebel et al., 2003;Roebroeck et al., 2005; IntroductIon Goal-oriented behavior requires outcome monitoring. Many studies in non-human primates described the role of the saliency/ reward pathway, involving the ventral tegmental area (VTA) and substantia nigra (SN), in outcome and error processing (Montague et al., 1996;Schultz et al., 1997;Holroyd and Coles, 2002;Schultz, 2002;Fiorillo et al., 2003;Bayer and Glimcher, 2005). Neurons in the VTA/SN increase activity to an unexpected reward and decrease activity to a missing reward. These roles of VTA/SN in error-related cognitive processes have recently been substantiated in humans (D'Ardenne et al., 2008;Carter et al., 2009;Duzel et al., 2009).
Anatomical and electrophysiological studies converged to suggest a function of the epithalamus/habenula in regulating outcome-related signals in the VTA/SN (Matsuda and Fujimura, 1992;Scheibel, 1997;Ji and Shepard, 2007;Hikosaka, 2007, 2009;Hikosaka et al., 2008;Morissette and Boye, 2008). For instance, neurons in the habenula and midbrain were each excited and inhibited by no-reward-predicting targets; and electrical stimulation of lateral habenula (LH) decreased activity in the dopaminergic neurons (Matsumoto and Hikosaka, 2007).
In an earlier functional magnetic resonance imaging (fMRI) study of humans performing a target prediction task, hemodynamic responses were observed in the anterior cingulate cortex, inferior anterior insula, and habenula during negative feedback (Ullsperger and von Cramon, 2003). The authors suggested a role Abler et al., 2006;Stilla et al., 2007;Deshpande et al., 2008;Duann et al., 2009;Sato et al., 2009;Ide and Li, 2011). We thus used GCA to ascertain the direction of connectivity between the habenula and VTA/SN. Furthermore, with mediation analysis we established that the habenula projected directly to the VTA/SN.

MaterIals and Methods subjects and behavIoral task
Fifty-nine healthy adult subjects (30 men, 22-45 years of age) participated in the study according to a protocol approved by Yale University Human Investigation Committee.
We employed a simple reaction time task in this stop signal paradigm (Logan et al., 1984;Li et al., 2006Li et al., , 2008aLi et al., , 2009). There were two trial types: "go" and "stop," randomly intermixed. A small dot appeared on the screen to engage attention at the beginning of a go trial. After a randomized time interval (fore-period) between 1 and 5 s, the dot turned into a circle (the "go" signal), prompting the subjects to quickly press a button. The circle vanished at a button press or after 1 s had elapsed, whichever came first, and the trial terminated. A premature button press prior to the appearance of the circle also terminated the trial. Approximately three quarters of all trials were go trials. The remaining one quarter were stop trials. In a stop trial, an additional "X," the "stop" signal, appeared after and replaced the go signal. The subjects were told to withhold button press when they saw the stop signal. Likewise, a trial terminated at button press or when 1 s had elapsed since the appearance of the stop signal. The stop signal delay (SSD) -the time interval between the go and stop signal -started at 200 ms and varied from one stop trial to the next according to a staircase procedure: if the subject succeeded in withholding the response, the SSD increased by 64 ms; conversely, if they failed, SSD decreased by 64 ms (Levitt, 1971). There was an inter-trial-interval of 2 s. Subjects were instructed to respond to the go signal quickly while keeping in mind that a stop signal could come up in a small number of trials. All participants had a practice session outside the scanner and completed four 10-min runs of the task in the scanner. Depending on the actual stimulus timing (trials varied in fore-period duration) and speed of response, the total number of trials varied slightly across subjects in an experiment. With the staircase procedure we anticipated that the subjects would succeed in withholding their response in approximately half of the stop trials.
As a control experiment, we also imaged 30 subjects during a 10-min resting state session, in which subjects were instructed to stay awake and relaxed, with their eyes closed (Duann et al., 2009).

IMagIng protocol
Conventional T 1 -weighted spin echo sagittal anatomical images were acquired for slice localization using a 3-T scanner (Siemens Trio). Anatomical images of the functional slice locations were next obtained with spin echo imaging in the axial plane parallel to the AC-PC line with TR = 300 ms, TE = 2.5 ms, bandwidth = 300 Hz/ pixel, flip angle = 60°, field of view = 220 mm × 220 mm, matrix = 256 × 256, 32 slices with slice thickness = 4 mm and no gap. Functional, BOLD signals were then acquired with a singleshot gradient echo echoplanar imaging (EPI) sequence. Thirty-two axial slices parallel to the AC-PC line covering the whole brain were acquired with TR = 2,000 ms, TE = 25 ms, bandwidth = 2,004 Hz/pixel, flip angle = 85°, field of view = 220 mm × 220 mm, matrix = 64 × 64, 32 slices with slice thickness = 4 mm and no gap. Three hundred images were acquired in each session. spatIal pre-processIng of braIn IMages Data were analyzed with Statistical Parametric Mapping 8 (Wellcome Department of Imaging Neuroscience, University College London, UK). Images from the first five TRs at the beginning of each trial were discarded to enable the signal to achieve steady-state equilibrium between RF pulsing and relaxation. Images of each individual subject were first corrected for slice timing, realigned (motion-corrected) and unwarped (Andersson et al., 2001;Hutton et al., 2002). A mean functional image volume was constructed for each subject for each run from the realigned image volumes. These mean images were co-registered with the high resolution structural image and then segmented for normalization to an Montreal Neurological Institute (MNI) EPI template with affine registration followed by non-linear transformation Friston, 1999, 2005). The normalization parameters determined for the mean functional volume were then applied to the corresponding functional image volumes for each subject. Finally, images were smoothed with a Gaussian kernel of 6 mm at full width at half maximum.

general lInear ModelIng
We followed our previous studies in the statistical modeling of imaging data (Li et al., 2006(Li et al., , 2008a. Briefly, four trial types were distinguished: go success (G), go error (F), SS, and SE trials. A statistical analytical design was constructed for each individual subject, using the general linear model (GLM) with the onsets of go signal in each of these trial types convolved with a canonical hemodynamic response function (HRF) and with the temporal derivative of the canonical HRF entered as regressors in the model (Friston et al., 1995). We entered reaction time (RT) and SSD as parametric modulators for go and stop trials, respectively, in the GLM. Realignment parameters in all six dimensions were also entered in the model. The data were high-pass filtered (128 s cutoff) to remove low-frequency signal drifts, and serial autocorrelation caused by aliased cardiovascular and respiratory signals was corrected by a first-degree autoregressive or AR (1) model (Friston et al., 2000;Della-Maggiore et al., 2002). Across subjects, there were: 281.6 ± 19.6 G trials, 11.2 ± 12.0 F trials, 47.2 ± 4.7 SS trials, and 40.8 ± 7.2 SE trials. We did not include F trials in the current analyses, because they comprised less than 3% of all trials.
In the first-level analysis, we constructed for each individual subject a contrast SE > SS in order to identify regional brain activations associated with error detection (Li et al., 2008b). Our previous work showed that the contrasts SS > G and SE > G activated brain regions that overlapped those of SE > SS, including the anterior cingulate cortex (both SS > G and SE > G) as well as the thalamus and midbrain (SE > G; Li et al., 2008b). Thus, to demonstrate the specificity of the contrast SE > SS in PPI, we examined SE > G and SS > G for comparison (see below). The contrast images (con) of the first-level analysis were used for random-effect analysis (Penny et al., 2004) to obtain group T maps using a one-sample t test. Brain regions were identified using an atlas (Duvernoy, 2003;Mai et al., 2008). time series were entered into multivariate autoregressive (MAR) modeling (Harrison et al., 2003). We used the Akaike Information Criterion (AIC), which imposes a complexity penalty on the number of parameters and avoids over-fitting of the data (Akaike, 1974). The application of MAR modeling required that each ROI time series was covariance stationary, which we examined with the Augmented Dickey Fuller (ADF) test (Hamilton, 1994). The ADF test verified that there was no unit root in the modeled time series. The residuals of MAR modeling were used to compute the Granger causality measures (F values) of each possible connection between ROIs. Since MAR modeling often involves highly interdependent residuals (Deshpande et al., 2009), we used permutation resampling (Hesterberg et al., 2005;Seth, 2010) to obtain an empirical null distribution of no causality, as suggested in Roebroeck et al. (2005), in order to estimate the F critical , and assess the statistical significance of Granger causality measures. With resampling, we produced surrogate data by randomly generating time series with the same mean, variance, autocorrelation function, and spectrum as the original data (Theiler et al., 1992), as implemented in previous EEG (Kaminski et al., 2001;Kus et al., 2004), and fMRI studies (Deshpande et al., 2009). We used binomial test to assess statistical significance in group analysis (Duann et al., 2009;Uddin et al., 2009); for each connection, we counted the number of subjects that had F > F critical (i.e., significant connection) and estimated its significance using a binomial distribution with parameters n = 59 trials and p = q = 0.5 (same probability to observe a connection or not). Multiple comparisons were corrected for false discovery rate (FDR; Genovese et al., 2002).
We applied the same multivariate GCA procedures to resting state data of the 30 subjects, following our previous work (Duann et al., 2009), as an additional control for false positive connectivities. The absence of functional connectivity in the resting data would suggest that the task-related connectivity is not an artifact of HRF variability across the brain (see next section for more details).

GCA: methodological considerations
Although some investigators argued the importance of causality based on temporal precedence and the utility of GCA in connectivity analyses (e.g., Roebroeck et al., 2009), others discussed the limitations of GCA (e.g., Friston, 2009). As detailed in a recent review of GCA in neuroimaging (Bressler and Seth, 2010), the utility of Granger causality measures depends on successfully estimating autoregressive (AR) models of stochastic processes. Successful applications of GCA to fMRI data (Roebroeck et al., 2005;Stilla et al., 2007;Bressler et al., 2008;Deshpande et al., 2009;Duann et al., 2009;Kayser et al., 2009;Ide and Li, 2011) appeared to have some elements in common. First, it is crucial that the modeled time series are wide-sense stationary (WSS; i.e., they have constant mean and variance). Otherwise, non-stationary time courses are known to produce spurious regression results (Granger and Newbold, 2001). Second, it is also important to have a number of observations (time points) adequate to estimate the AR model coefficients. In the current GCA modeling, by concatenating BOLD time series across four sessions for each individual (a total of 1,180 time points) after de-trending and normalization, we obtained time series (averaged inside each ROI) that were sufficiently long and covariance stationary, the latter verified by the ADF test. Third, we applied spatial psychophysIologIcal InteractIon We used PPI to describe how functional connectivity between brain regions was altered as a result of psychological context (Friston et al., 1997). We hypothesized that the habenula showed greater PPI with the VTA and SN during SE compared to SS trials. With the contrast SE > SS, we identified a mask of the habenula comprising two symmetric spheres each of 6 mm in radius and centered at MNI coordinates [−1, −25, 1] and [1, −25, 1]. This mask was well within the area identified as habenula in Ullsperger and von Cramon, 2003 (Talairach and approximate MNI coordinates: [−5/6, −25, 8] and [−6/6, −25, 5]). On the basis of the GLM, we extracted the time series of the first eigenvariate of the BOLD signal of the habenula for each individual subject. The eigenvariate value, inside a region of interest (ROI), corresponds to the average BOLD signal weighted by the voxel significance, and it is more robust to outliers (Gitelman et al., 2003). This time series constituted the physiological variable. The time series were de-convolved to remove the effects of HRF, multiplied by the psychological variable (SE > SS, i.e., "1" for SE and "−1" for SS conditions), and re-convolved with the canonical HRF to obtain the interaction term or PPI variable (Gitelman et al., 2003). The three variables were entered as regressors in a whole-brain GLM. PPI analysis was performed for each individual subject, and the resulting positive contrast images (i.e., "1" for the PPI regressors) were used in random-effect group analysis (Penny et al., 2004). Group results were reported for p < 0.001, uncorrected. To check whether PPI results were specific to the contrast SE > SS, we also used SE > G and SS > G as psychological variables of PPI for comparison.
To ascertain the specificity of brain regions of PPI with the habenula during error processing, we performed another PPI analysis with a "control" seed region, involving the error activated thalamic cluster (MNI coordinates [14,5] and [−14, −20, 10], 5,600 mm 3 ), exclusively masked by a larger area encompassing the habenula (two spheres of 10 mm in radius and centered at [−1, −25, 1] and [1, −25, 1]). This mask was to ensure the exclusion of the habenula from the control seed region.
For the regression slope analysis, multivariate GCA, and mediation analysis, we referred to the habenula and the regions identified with PPI as the ROIs or only the latter brain regions as ROIs.

MultIvarIate granger causalIty analysIs
The analysis of PPI identified areal interaction but did not specify the direction of influence. In order to confirm our hypothesis that the habenula signals the VTA in error detection and not the other way around, we used a multivariate GCA (Stilla et al., 2007;Deshpande et al., 2009) to examine the direction of the influence between the ROIs of PPI and the habenula.
Multivariate GCA was implemented as in our previous work (Duann et al., 2009;Ide and Li, 2011). We considered two models; the first one included the habenula and all four regions identified from PPI as ROIs, and the second model included the habenula and three PPI regions (all except the globus pallidus; see below). Multivariate GCA was performed for individual subjects. For each subject and each ROI, a summary time series was computed by averaging across voxels inside the ROI for each time point. These average time series were concatenated across four sessions, after de-trending and normalization (Ding et al., 2000). Afterward, the pre-processed but not temporal filtering to the original BOLD signals because temporal (e.g., bandpass) filtering is known to introduce severe confounds in GCA of neuroimaging time series (Florin et al., 2010;Seth, 2010). These procedures were also successfully applied in our previous studies (Duann et al., 2009;Ide and Li, 2011).
An additional consideration is the effects of HRF variability and down-sampling of BOLD signals on autoregressive modeling. A popular approach is to use the "difference of influence" (DOI) between two regions (Roebroeck et al., 2005;Stilla et al., 2007) to ameliorate these effects. However, this approach is no longer valid for multivariate GCA. In such cases, a useful practice is to analyze Granger causality during different experimental conditions (e.g., by studying both task and resting conditions; Duann et al., 2009;Kayser et al., 2009), since the effects of HRF variability and signal down-sampling are not expected to vary across conditions.
A few studies in the literature presented less successful results from GCA. In a study of simultaneous electroencephalographic recordings and fMRI in rats, David et al. (2008) showed spurious connectivities derived from GCA due to HRF variability. However, one should note that the HRF variability was outside normal physiological range. Furthermore, in a recent investigation using simulated BOLD signals by convolving a standard HRF with local field potentials recorded from macaque cortex, Deshpande et al. (2010) showed that, even considering real and normal range of HRF variability (Handwerker et al., 2004) and a signal-to-noise ratio of 1 unit and a TR = 2 s, GCA reliably detected neuronal delays around 700 ms (Deshpande et al., 2010). Witt and Meyerand (2009) reported poor performance of GCA but it is not clear whether these experiments were biased because of simulated fMRI time series generated using Dynamic Causal Modeling (Friston, 2009) or whether these simulated time series were covariance stationary.
Taken together, although GCA presents some technical challenges, we believe that, when it is carefully applied and done without over-interpretation of the results, GCA is a useful exploratory tool to delineate effective connectivity of the complex human brain.

MedIatIon analysIs
We performed mediation analyses to further characterize functional connectivity between the ROIs (MacKinnon et al., 2007), using the toolbox M3, developed by Tor Wager and Martin A. Lindquist 1 . Mediation analyses are widely used in social and economic research to examine whether a relationship between two variables is mediated by an intervening variable (Maccorquodale and Meehl, 1948;Baron and Kenny, 1986). It was successfully applied to fMRI of emotion regulation (Wager et al., 2008;Lebrecht and Badre, 2008) and more recently to analyses of functional connectivity (Hare et al., 2010). In a mediation analysis, the relation between the independent variable X and dependent variable Y, i.e., X → Y, is tested to see if it is significantly mediated by a variable M. The mediation test is performed by employing three regression equations (MacKinnon et al., 2007): where a represents X → M, b represents M → Y (controlling for X), c′ represents X → Y (controlling for M), and c represents X → Y. In the literature, a, b, c, and c′ were referred as path coefficients or simply paths (MacKinnon et al., 2007;Wager et al., 2008), and we followed this notation. Variable M is said to be a mediator of connection X → Y, if (c − c′) is significantly different from 0, which is mathematically equivalent to the product of the paths a × b (MacKinnon et al., 2007). If the product a × b and the paths a and b are significant, one concludes that X → Y is mediated by M. In addition, if path c′ is not significant, it indicates that there is no direct connection from X to Y and that X → Y is completely mediated by M. Note that path b is the relation between Y and M, controlling for X, and it should not be confused with the correlation coefficient between Y and M.

Mediation analyses: methodological considerations
As with other methods based on structural equation models, one assumed that all relevant variables are included in the analysis; i.e., one could not rule out the existence of mediating factors not tested in the model (Lebrecht and Badre, 2008). In addition, mediation analysis is only valid upon correct specification of the causal orders (MacKinnon et al., 2007). We believe that these limitations were addressed to a significant extent in the current study: wholebrain PPI identified all relevant ROIs functionally connected to the habenula during errors; and multivariate GCA provided important information regarding causal orders. Finally, as pointed out by Wager et al. (2008), an additional limitation of using mediation analysis in fMRI is that models are made on the basis of naturally occurring variance over subjects, and thus conclusions are made with the assumption that inter-subject variability does not affect the coupling between dependent variables. This restriction also applies to the study. One could control the variability by simply removing the regression outliers (Chatterjee and Hadi, 1986) or, alternatively, develop multilevel mediation models that consider the mediation path coefficients as random effects (MacKinnon et al., 2007).

results braIn regIons of ppI wIth the habenula
With general linear modeling we examined and confirmed errorrelated regional brain activations during the stop signal task (Li et al., 2008b). Compared to SS, SE trials evoked greater activations in the medial frontal cortex including the dorsal anterior cingulate cortex (dACC) and pre-supplementary motor area (preSMA), as well as the anterior inferior insulas, thalamus, habenula, and structures in the midbrain, at p < 0.05, corrected for family-wise error (FWE) of multiple comparisons (Figure 1; Table 1).
In PPI, we identified brain regions that were functionally connected with the habenula. Compared to SS trials, SE trials evoked greater PPI with the habenula in the VTA/SN, bilateral anterior inferior insula, amygdala, and the internal segment of the globus pallidus (GPi), at p < 0.001, uncorrected (Figure 2; Table 2). To test our hypothesis targeting the VTA/SN, we also performed a ROI analysis, using small volume correction for a spherical mask centered at [2, −16, −10] and 6 mm in radius (Carter et al., 2009). The results identified a peak of activation at [6, −16, −14], p < 0.005, corrected for FWE. In addition, on the basis of the literature (Seeley et al., 2007;Hikosaka et al., 2008;Haber and Knutson, 2010), Because voxel smoothing in image pre-processing diminished spatial specificity of the results, we examined brain regions of PPI with the thalamic cluster "minus the habenula" as a control. The results showed that, at the same statistical threshold, regions showing PPI with the thalamus included the occipital cortices, fusiform gyrus, temporal cortex, insula, supplementary motor area, and prefrontal cortex (Figure 5; Table 3). No activation was observed in the VTA/SN region (p < 0.01, uncorrected). Thus, the PPI of the VTA/SN, bilateral amygdala, and the GPi appeared to be specific to the habenula.

granger causalIty analysIs
With multivariate GCA, we determined the direction of influence between the habenula and ROIs of PPI (Figure 6). In one model we included the habenula, VTA/SN, GPi, bilateral amygdala, and insula. To minimize model complexity and facilitate interpretation of the results, we combined bilateral amygdala as a single ROI. Similarly, we combined bilateral insula as a single region. For each individual subject we tested all possible connections between regions at p < 0.05, corrected for FDR, against an empirical statistical null distribution (see Materials and Methods). In a binomial test (p < 0.05) for the group, the results showed significant projections from the habenula to the VTA/SN, insula, and amygdala, but not to the GPi ( Figure 6A). Thus, in a second model, we excluded the GPi, and again obtained significant projections from the habenula to the VTA/SN but not from the VTA/SN to the habenula ( Figure 6B). Thirty-eight of the 59 participants also showed a significant projection from the amygdala to the VTA, compared to 44 with projection from the habenula to the VTA/SN. In contrast to these task-related connectivity patterns, no significant Granger causality was observed (binomial test, p > 0.50, in both models) for any of the connections in the resting state time series, providing evidence that these task-related connections were not spurious, such as resulting from HRF variability across the brain. the insula, amygdala, and GPi were also significant at p < 0.05, corrected for FWE, after small volume correction using anatomical masks from the Automated Anatomical Labeling atlas (Tzourio-Mazoyer et al., 2002). Conversely, no brain regions showed greater connectivity with the habenula during SS compared to SE trials (p < 0.01, uncorrected).
In comparison, PPI with SE > G and SS > G as psychological variables each revealed activations in the bilateral insula (SE > G) and precuneus (G > SE), and the right insula (SS > G) as well as posterior cingulate cortex and the precuneus (G > SS), p < 0.001, uncorrected (Figures 3 and 4). No foci were observed in the area of the VTA/SN or the amygdala (p < 0.01, uncorrected).

MedIatIon analysIs
We performed mediation analyses to further characterize the functional connectivity between the habenula, amygdala, and VTA/SN during error processing. In particular, we tested the hypothesis that the projection amygdala → VTA/SN was mediated by the habenula. We derived for each individual subject the effect-size of SE > SS of the three ROIs for a single-level  Brett et al., 2002 2 ). We tested two models. In the first model, the effect-sizes of the amygdala and VTA/SN, habenula were set as X, Y, and mediator variable M, respectively. In the second model, we tested if the projection from habenula (X) to VTA/SN (Y) was mediated by amygdala (M). Figure 7 and Table 4 summarized the results. The results of the first model confirmed the hypothesis of the habenula mediating the connection from amygdala to VTA/SN during error processing. The results of the second model confirmed the lack of a projection from the amygdala to the VTA/SN.

dIscussIon
The VTA/SN, bilateral amygdala, insula, and GPi showed greater PPI with the habenula during SE as compared to SS trials. In contrast, except for the left insula, none of these areas showed a PPI with the thalamic cluster. Thus, with the exception of the left insula, these brain regions altered activation to SE as compared to SS trials in specific association with the habenula in the stop signal task. Furthermore, the connectivity between the habenula and VTA/SN was contrast-specific as PPI of SS > G or SE > G did not reveal this connection.
As described in the above, PPI analyses did not provide information about the direction of influence between brain regions. Thus, we performed GCA, which confirmed a feedforward connection from the habenula to VTA/SN but not vice versa. The amygdala also showed a significant projection to the VTA/SN in addition to a bidirectional connection with the habenula, broadly consistent Unlike single-unit recordings, which could determine the "sign" of influence of neuronal activities in one brain region on those in another, the BOLD signals in fMRI represented an indirect measure of neural activity (Logothetis, 2008). Thus, although the current findings demonstrated a feedforward connectivity between the habenula and the VTA/SN, they did not indicate whether the influence is excitatory or inhibitory. Indeed, Ji and Shepard (2007) and Hikosaka et al. (2008) showed that neurons in the LH might exert inhibitory effects on neurons in the VTA/SN by exciting local GABAergic neurons. Furthermore, dopaminergic neurons in the VTA/SN may respond to stimuli with both positive and negative motivational value (Matsumoto and Hikosaka, 2009). These issues may lead one to consider whether the functional connectivity between the habenula and the midbrain had more to do with saliency than the motivational valence of errors. On the other with evidence of anatomical projections from the amygdala to the VTA (Kaufling et al., 2009) and the role of amygdala in reward and saliency processing (Baxter and Murray, 2002;Etkin et al., 2006;Murray, 2007;Haber and Knutson, 2010;Linke et al., 2010; see also Delgado et al., 2008 for a review). On the other hand, Granger causality did not describe event-related relationship between time series. Thus, we used mediation analyses to further characterize the functional connectivities between the habenula, amygdala, and VTA/SN during error processing. The results supported the hypothesis that the error-related connectivity between the amygdala and VTA/SN was largely mediated by the habenula and that the error-related connectivity between the habenula and VTA/SN was not mediated by the amygdala. Taken together, the results of these complementary analyses support a robust error-related signal from the habenula to VTA/SN during the stop signal task.  Table 3. Note that, except for the left insula, the brain regions of PPI with the thalamus cluster did not overlap those of PPI with the habenula (Figure 2).

Methodological considerations
Some methodological issues need to be considered. First, the volume of the VTA is approximately 60 mm 3 (Paxinos and Huang, 1995); and the volume of the SN, pars compacta, is approximately 190 mm 3 , estimated based on atlas plates (Mai et al., 2008 (Wittmann et al., 2005). The same localization issue applies to the habenula, a small structure with a size of approximately 3 mm × 3 mm × 6 mm (Salas et al., 2010). We sought to obtain the functional time course reflecting the habenula activity, by computing the average BOLD signal of a mask with a volume of approximately 1,344 mm 3 or ∼29 voxels. Note that this mask included voxels in the white matter and CSF and raised the question whether the functional connectivities were specific to the habenula. We wish to address this issue from two perspectives. First, as described in the Section Materials and Methods, aliased cardiovascular and respiratory signals were removed with high pass filtering and AR (1) in analyses (Friston et al., 2000;Della-Maggiore et al., 2002). Second, the PPI of the habenula mask showed results that were anatomically specific while the control mask (which also contained white matter and CSF) did not. We feel that signal artifacts of the CSF or white matter are unlikely to account for the current results. In addition, we computed PPI for a "control" region that excluded this habenula mask for comparison.
Second, GCA per se does not involve modulation of the connectivity by experimental conditions, as does PPI analysis. Thus, GCA alone does not guarantee "causality" between regional brain activations in response to specific events in the cognitive task. That is, it was possible that the causalities observed between the habenula, amygdala, and the VTA/SN was due to signals unrelated to error processing. For this reason, we used mediation analyses to examine the correlation of activity between these brain regions. In particular, compared to the habenula, the amygdala hand, both SE and SS trials are more salient than G trials. Our supplementary results showing a largely lack of activation in the VTA/SN during PPI with SE > G and SS > G suggested that the feedforward influence of the habenula on VTA/SN was contingent upon a contrast between negative and positive motivational value.
In rodent studies, Jhou et al. (2009) reported the mesopontine rostromedial tegmental (RMTg) nucleus mediating the projections of the LH to VTA/SN. In the PPI analyses, we did not observe any significant clusters in the brain stem that might correspond in location to the RMTg nucleus. In contrast, we observed activation in a cluster around the pretectal area (Figure 2, z = −14), consistent with an earlier work showing projection from the superior colliculus and pretectal area to SN during the detection of salient visual events (Comoli et al., 2003). Although the GPi was functionally connected to the habenula during PPI, neither GCA or regression slope analysis could further confirm its role in this circuit. The latter result thus needs to be reconciled with findings from single-unit recordings of monkeys showing habenula-projecting neurons in the GPi (Hong and Hikosaka, 2008).  results at (B) show that the habenula, insula, and amygdala are bidirectionally connected, while the habenula and amygdala unilaterally project to the VTA/SN. GCA was performed and evaluated for each connection of each individual subject, at p < 0.05, corrected for FDR. Group results were obtained using a binomial test for each connection (p < 0.05). The numbers next to the arrows indicate the number of subjects out of 59 who have that connection.  humans. These new findings extend recent work in rodents and non-human primates on the role of this subcortical circuit in error detection.

acknowledgMents
This study was supported by NIH grants R01DA023248 (Chiang-Shan R. Li), K02DA026990 (Chiang-Shan R. Li), R21AA018004 (Chiang-Shan R. Li), and a Clinical Translational Science Award (NIH-UL1RR024139, Robert Sherwin). The NIH had no further role in study design; in the collection, analysis, and interpretation of data; in the writing of the report; or in the decision to submit the paper for publication. We thank Drs. Wager and Lindquist for making the toolbox of mediation analysis freely available on the web. We also thank Drs. Sheng Zhang, Sien Hu, Olivia Hendrick, and Sarah Bednarski for their many helpful discussions.
showed a less significant projection (Granger causality) to the VTA/ SN, raising the possibility that the habenula may mediate the projection from the amygdala to VTA/SN during error processing. The results of mediation analyses substantiated this hypothesis. Error-related activity of the habenula was significantly correlated with the VTA/SN even with the amygdala as an intervening variable; in contrast, amygdala correlated with the VTA/SN largely through its connection with the habenula. Finally, a contrast of SE > SS involved differences in motor responses in addition to error processing ( Table 1). Although to our knowledge the habenula has never been implicated in motor control, the current results need to be considered in the context of this limitation.

conclusIon
The complementary results of psychophysiological, regression slope, Granger causality, and mediation analyses delineated a feedforward connection from the habenula to the VTA/SN in Frontiers in Human Neuroscience www.frontiersin.org