|Home | About | Journals | Submit | Contact Us | Français|
Granger causality analysis of functional magnetic resonance imaging (fMRI) blood-oxygen-level-dependent signal data allows one to infer the direction and magnitude of influence that brain regions exert on one another. We employed a method for upsampling the time resolution of fMRI data that does not require additional interpolation beyond the interpolation that is regularly used for slice-timing correction. The mathematics for this new method are provided, and simulations demonstrate its viability. Using fMRI, 17 snake phobics and 19 healthy controls viewed snake, disgust, and neutral fish video clips preceded by anticipatory cues. Multivariate Granger causality models at the native 2-sec resolution and at the upsampled 400-ms resolution assessed directional associations of fMRI data among 13 anatomical regions of interest identified in prior research on anxiety and emotion. Superior sensitivity was observed for the 400-ms model, both for connectivity within each group and for group differences in connectivity. Context-dependent analyses for the 400-ms multivariate Granger causality model revealed the specific trial types showing group differences in connectivity. This is the first demonstration of effective connectivity of fMRI data using a method for achieving 400-ms resolution without sacrificing accuracy available at 2-sec resolution.
In the two decades since the blood-oxygen-level dependent (BOLD) effect for functional magnetic resonance imaging (fMRI) was first described, there has been a rapid proliferation of improved methods for data acquisition and statistically rigorous analysis. This is particularly true of functional connectivity methods for characterizing the interactions between areas distributed throughout the brain (Friston et al., 1995, 1997; Gitelman et al., 2003). Effective connectivity measures go a step further by directly addressing causation (Deshpande et al., 2009; Friston et al., 2003; Goebel et al., 2003).
Since cause should precede effect, the key aspect to inferring causation is temporal lag. Granger causality utilizes lagged time series of one or more variables to determine those that are useful in predicting future values of one or more variables (Granger, 1969). A problem with bivariate Granger models is that a third variable may be responsible for the presumed causal influence identified between two variables. This study employed a multivariate approach to guard against the risk of finding false causation, although it is still susceptible to variables not included in the model.
Another obstacle for inferring direct causation from fMRI data is that multisynaptic events can occur in the time needed to collect a single-brain volume. Typically, time of repetition (TR) for acquisition of a whole-brain functional volume is in the range of 2sec, whereas synaptic events can occur on the scale of hundreds of milliseconds or faster. We devised a method to increase the temporal resolution of fMRI data that capitalizes on standard acquisition sequences using interleaved slices. Commonly implemented at an early stage of fMRI data preprocessing, slice-timing correction adjusts for the time offset between these interleaved slices by interpolating between consecutive acquisitions of an individual slice to estimate each voxel's value at a particular instant in time. A preprocessing step that is currently included in all major fMRI processing packages (Sladky et al., 2011) adjusts each interleaved slice to a single time point in the volume, such that all other slices are interpolated to that time point. The method proposed here instead adjusts each slice to multiple time points at equivalent increments within a single volume, thus upsampling the data during slice-timing correction. Upsampling later in data processing requires interpolating data to a single time point during slice-timing correction, and then interpolating data a second time to achieve a finer resolution. The first method should retain information that is lost by the second method and, therefore, capture more of the dynamics in the original time series. Also, the proportion of data interpolated using this new upsampling method is identical to the proportion of data interpolated for conventional fMRI analysis (e.g., 97% for 30-slice acquisition).1
Granger causality is an effective and proven methodology for increasing the information and insight provided by fMRI (Bressler and Seth, 2010; Rogers et al., 2010). Recent reports comparing Granger causality and dynamic causal modeling have identified mutual strengths as well as advantages and disadvantages of each (David et al., 2008; David, 2009; Friston 2009a,b; Roebroeck et al., 2009a,b; Zou et al., 2009). Granger causality is more effective at exploration and model development, whereas dynamic causal modeling is particularly well suited for testing models (Friston, 2009b; Roebroeck et al., 2009a,b). In addition, Granger causality is superior with long data sets, such as in the current study (Zou et al., 2009). Multiple putative disadvantages of Granger causality (Friston, 2009a) have been retracted (Friston, 2009b) or resolved with advances in the application of Granger causality to fMRI analysis (Deshpande et al., 2010a; Havlicek et al., 2010; Liao et al., 2009; Marinazzo et al., 2010; Sato et al., 2010; Sladky et al., 2011; Zhou et al., 2009a,b).
In the present study, we directly compared multivariate Granger causality analysis models at the native 2-sec resolution with models at an upsampled 400-ms resolution. A context-dependent multivariate Granger causality analysis was also performed. Subjects with snake phobia and healthy volunteers underwent an event-related fMRI paradigm that involved viewing emotional video clips preceded by anticipatory cues: S cues always preceded videos of moving snakes, F cues always preceded videos of swimming fish, and D cues always preceded videos of disgusting scenes. All analyses were conducted on anatomically defined regions of interest (ROIs) implicated in prior research on anxiety and emotion, including the amygdala, hippocampus, anterior insula, and anterior cingulate cortex (ACC) (Dalgleish, 2004; Davis et al., 2010; Etkin and Wager, 2007; Nitschke et al., 2006, 2009; Phillips et al., 2003a,b).
The rationale for the development of this method is that improving temporal resolution might retain information that is otherwise lost at native resolution and, thus, better reflect the dynamics of the original data, which was investigated for both simulated data and real-world data. We further hypothesized that additional information at 400-ms resolution would increase the number of significant influences detected between brain regions.
The 19 snake phobic participants (15 females, mean age 21.8, range 18–46) and 17 healthy controls (7 females, mean age 28.6, range 19–58) were right-handed and neurologically healthy. Snake phobics met criteria for DSM-IV diagnosis of specific phobia of snakes and had no other current or past diagnoses, as assessed by the Structured Clinical Interview for the DSM-IV (SCID; First et al., 2002). Controls were healthy volunteers with no history of mental disorder, as assessed by the SCID.2
Whole-brain anatomical and functional images were acquired on a GE 3.0T scanner and processed using Analysis of Functional NeuroImages (AFNI) (Cox, 1996) using standard procedures described in Supplementary Methods (Supplementary Data are available online at www.liebertpub.com/brain) and prior publications (Nitschke et al., 2009; Sarinopoulos et al., 2010). Functional scans were acquired as 30 interleaved sagittal T2*-weighted echo planar image slices (TR/time of echo=2000/30ms, flip angle=90°, number of excitations=1, field of view=240mm, matrix=64×64, slice thickness/gap=4.5/0.5mm, in-plane resolution=3.75×3.75mm). The functional paradigm was an emotional anticipation task involving simple cues predicting video clips (Fig. S1; Supplementary Methods). In the controllable condition (but not the uncontrollable condition), participants could prevent the video presentation if their responses to a target were fast enough. Based on prior research on anxiety and emotion (Dalgleish, 2004; Davis et al., 2010; Etkin and Wager, 2007; Nitschke et al., 2006, 2009; Phillips et al., 2003a,b), we selected a set of 13 anatomical ROIs: bilateral dorsal amygdala, bilateral ventral amygdala, bilateral hippocampus, bilateral anterior insula, bilateral subgenual ACC, bilateral supragenual ACC, and pregenual ACC (Fig. S2; Supplementary Methods).
Fourier interpolation at the voxel level with AFNI slice-timing correction (3dTshift) was utilized to upsample data without any additional interpolation than that which is regularly used in fMRI preprocessing. Rather than correcting to a single time point per volume as is typically done (e.g., first slice of a volume or halfway point of a volume), data were corrected to slices at 400-ms increments. Since data were acquired as 30 interleaved slices, the 400-ms increments were demarcated by slices 0, 12, 24, 7, and 19. The slice nnn option of 3dTshift was used to have the slice timing corrected to these slices (i.e., nnn=0, 12, 24, 7, and 19). The 400-ms increments resulted in five separate whole-brain volumes per TR that were then run through the rest of preprocessing (motion correction via volume registration, fieldmap correction, and conversion to percent signal change [PSC]). Of note, an ongoing question in fMRI analysis is whether slice-timing correction should be performed before or after motion correction (Roche, 2011). We implemented slice-timing correction before motion correction so that the Fourier interpolation equation generated for upsampling would be based on original data rather than based on data interpolated by motion correction. As such, motion correction was performed at the finer temporal resolution. Applying motion correction after our upsampling slice-timing procedure resulted in the upsampled data being assigned to the correct voxel. 3dmaskave was then used to extract mean PSC time series for each ROI from each of the five volumes. The output of 3dmaskave was five time series at 2-sec resolution that were 400ms apart. These five time series were then assembled into a single time series at 400-ms resolution. Equations for upsampling and slice timing correction are provided in Supplementary Methods.
Resolutions other than 400ms were considered for upsampling via the slice-timing procedure described above. Acquiring 30 slices per whole-brain volume allows resolutions of 67ms, 133ms, 200ms, 333ms, 400ms, 667ms, and 1000ms, as determined by dividing the TR (2000ms) by all integers divisible into 30 (1, 2, 3, 5, 6, 10, 15, 30). A cost to increasing resolution is the corresponding increase in the number of comparisons being made for an equivalent time frame, which results in increased stringency of the Bonferroni correction and the reduction of statistical power. A 400-ms model was selected over the models with even finer resolution for two reasons: (1) these concerns about statistical power, and (2) computational limits, which prevented the models at finer resolution from covering a time span corresponding to at least two lags (4sec) of the standard model at 2-sec resolution to allow direct comparisons between the models.3
For the 400-ms model, 10 lags at 400-ms increments (400–4000ms) were modeled simultaneously by vector auto regression.
where q is the number of ROIs, p is the number of lags, cg,q,h is the path coefficient from region g to region q at lag h. The time series was sampled at 400ms (7235 time points), which is used by 1dGC.R in AFNI (Chen et al., 2009; http://afni.nimh.nih.gov/sscc/gangc/1dGC) to solve for the influences that occur at 400, 800, 1200, 1600, 2000, 2400, 2800, 3200, 3600, and 4000ms. For the 2-sec model, three lags were used (2, 4, and 6sec). The number of lags for both models was determined using Akaike Information Criterion (AIC), with the 10 lags for the 400-ms model encompassing 4sec, which corresponds to the first two lags of the 2-sec model.4 Both models included the 13 ROIs and linear baseline terms. The mean of PSC time series values from each ROI was used in single-subject analyses. The path coefficients and t statistics from 1dGC.R at the single-subject level were then used by 1dGC.R in a mixed effects model to determine which paths were significant at the group level. 1dGC.R carries the correlation signs to the group level, thus indicating whether paths are positive or negative influences.5
The upsampling procedure described in “Data acquisition, processing, and upsampling” is a deterministic transformation (i.e., for a given input, this transformation will always produce the same output) of the original data and, thus, cannot artificially inflate significance. As long as a transformation is deterministic and is applied to all data equally, it is allowable statistically (assuming appropriate application of the subsequent statistic tests). As a two-stage analysis, everything done for first-level analysis (upsampling, Granger causality, even first-level statistics) can be considered a data transformation in preparation for second-level statistic analysis.
In addition to the above methods for context-independent Granger causality across the entire time series to test for significant pathways irrespective of trial type, context-dependent Granger causality was implemented to test effective connectivity for each of the specific trial types separately. Since 1dGC.R allows breaks in the input data, we separated our 400-ms time series for the entire experimental session into smaller time series for each of the six trial types (i.e., contexts): controllable snake, controllable fish, controllable disgust, uncontrollable snake, uncontrollable fish, and uncontrollable disgust. To achieve this, an in-house program (www.brainimaging.wisc.edu/~mcfarlin) written by D.R.M. used the event timing files for each trial type to generate two files. The first file listed the trial type for each time point, which was then combined with the time series files for each region using 1dcat. The time series was then broken into separate files for each of the six trial types (from the start of a cue until the next cue) with sed, using the information provided by the first file. The second file contained the duration of each individual trial, thereby denoting when breaks in the time series for each trial type occurred, which is critical for 1dGC.R modeling so that activity in one region at the end of one trial is not used to predict activity in another region at the start of the next trial. These steps resulted in six time series (one for each trial type) and information about trial durations that were then input into 1dGC.R to assess Granger causality separately for each of the six trial types (i.e., contexts).
These methods for context-dependent Granger causality were applied in two ways. First, to determine which trial types if any were driving the differences found in the context-independent analyses described above, context-dependent analyses were conducted in a post hoc fashion for only those pathways that were significant for the context-independent analyses. Second, context-dependent Granger causality was applied in an omnibus analysis across all regions and interactions for the 400-ms model. For data at native resolution, context-dependent Granger causality analysis was not possible because the number of data points per trial at a 2-sec resolution was too low for a model to converge.6
Second-level analyses were performed as one-sample t tests for each group, and as two-sample t tests for comparisons between groups. For the 400-ms models, the corrected threshold corresponding to p<0.05 using Bonferroni across the 13 ROIs was p<0.00003 (0.05/(13*13*10)). For the 2-sec models, the less stringent threshold of p<0.0001 (0.05/(13*13*3)) was applied, due to the fewer number of pathways tested across 3 lags as compared to the 10 lags of the 400-ms models. For the context-dependent Granger causality omnibus analysis applied to all regions and interactions for the 400-ms model, the appropriate Bonferroni correction was again p<0.00003. The false discovery rate (FDR) based on Bonferroni corrections was also reported for all analyses.7
Data were simulated with random Gaussian drift at 15Hz, corresponding to 67-ms resolution, and otherwise in accordance with the parameters described above for the real-world data. Superimposed on this random 15-Hz drift was a 10-lag, 2.5-Hz interaction network between 13 nodes corresponding to our 400-ms model with 13 ROIs (lags at 400, 800, 1200, 1600, 2000, 2400, 2800, 3200, 3600, and 4000ms). Activity for each node was calculated from the 10 previous lags at the 13 nodes multiplied by path coefficients for each lag and summed (Equation 3). To simulate fMRI acquisition of BOLD signal, each 15-Hz time series was convolved with a standard HDR and sampled from interleaved slices at 0.5Hz, which corresponds to the 2-sec TR used for the real-world data. This simulation used five interleaved slices to match the modal number of slices for the 13 ROIs described above; adjacent slices were separated by 1sec (i.e., slices 1, 3, and 5 were separated by 67ms, and slices 2 and 4 occurred 1sec after slices 1 and 3, respectively). The 0.5-Hz time series for each slice was upsampled with Fourier interpolation to 15Hz and phase-shifted to correct differences in time of acquisition, which corresponds to the standard slice-timing correction procedure implemented for the real-world data. Three simulated data sets were created:
Granger causality modeling was then implemented on these simulated data to test the ability to detect path coefficients at 400-ms resolution.
Upsampling at the time of slice timing correction (voxel-level interpolation) was closer to simulated data than upsampling after calculating the ROI average (ROI-level interpolation), t=2.1, p<0.028 (Figs. 1A & S7). Estimates for the voxel-level interpolation were 25% more accurate than ROI-level interpolation, demonstrating the benefit of interpolating the data for upsampling and slice-time correction simultaneously rather than interpolating the data a second time after slice-time interpolation has been completed. Implementing Granger causality on the simulated 400-ms data indicated that the output coefficients and corresponding t values correlated with the input path coefficients even at a single-subject level, r=0.24, p<0.0001, and r=0.59, p<0.0001, n=1690, respectively. This indicates that the strength of the Granger causality weights and t values increased as the ground truth connectivity of the simulated data increased.
For visualization purposes, Figure 1B illustrates that the time series are more dynamic for voxel-level interpolation using slice-timing correction to 400-ms increments than for ROI-level interpolation implemented after conventional slice-timing correction to the initial slice of each volume. The increased frequency of changes in the direction of the signal intensity provides evidence that more information about regional brain activity is available in the fMRI data than can be detected at native 2-sec resolution. Since Granger causality relies on the dynamics of a time series for finding predictive variables, providing the model with more dynamic data should increase the detection of predictive variables.
Supporting our hypothesis that finer temporal resolution would also increase sensitivity, the 400-ms model consistently outperformed the 2-sec model (Table 1). For Bonferroni-corrected data, the total number of influences found for the 400-ms model was 241, compared to 85 for the 2-sec model (and an intermediate 219 for an exploratory 1-sec model). Similarly, for FDR-corrected data, the total number of influences found for the 400-ms model was 693, compared to 179 for the 2-sec model (and an intermediate 615 for an exploratory 1-sec model). For the control group, the 400-ms model resulted in 136 influences (directional effective connections) among the 13 ROIs, as compared to 37 influences for the 2-sec model (χ2=57, df=1, p<0.000001). The pattern was the same for the phobic group, with 97 influences at 400ms and 47 at 2sec (χ2=17, df=1, p<0.00003). The assessment of group differences indicated eight influences for the 400-ms model and one for the 2-sec model (χ2=5.4, df=1, p<0.02). The increased number of significant connections in the 400-ms model does not result from the greater number of lags tested in that model as compared to the 2-sec model, because the Bonferroni correction procedure used accounts for the number of pathways tested (p<0.00003 for the 400-ms model vs. p<0.0001 for the 2-sec model; see Methods). The increased sensitivity of the 400-ms model was also observed when FDR correction was applied (control group, 393 for 400-ms model versus 83 for 2-sec model, χ2=202, df=1, p<0.000001; phobic group, 289 versus 95, χ2=98, df=1, p<0.000001; group differences, 11 versus 1, χ2=8.3, df=1, p<0.004).
Of 14 interregional connections in the control group for the 2-sec model using Bonferroni correction, 6 were between homologous regions, including right to left influences for the hippocampus and left to right influences for the anterior insula. Of 8 nonhomologous connections, 3 emanated from the left hippocampus and 3 from the left anterior insula (Fig. 2).
Of 25 interregional connections in the phobic group, 11 were between homologous regions, including the same hippocampus and anterior insula influences noted above for controls. Of 14 nonhomologous connections, 2 were from the right dorsal amygdala to the bilateral ventral amygdalae, and the remaining 12 involved the hippocampus, including two that were observed for controls (left hippocampus to left anterior insula and left anterior insula to right hippocampus).
The only significant difference between the groups for the 2-sec resolution model was from right hippocampus to right supragenual ACC at the 4-sec lag. The significant inhibitory influence of this pathway noted above in phobic subjects was not present in control subjects.
The 400-ms model revealed many rapid influences that were not observed for the 2-sec model (Fig. 3). The control group had 66 interregional connections for the 400-ms model using Bonferroni correction, and 32 occurred before the first lag of the 2-sec model. The phobic group had 26 interregional influences, 8 of which occurred before the first lag of the 2-sec model. Of particular note for both groups are bilateral influences of the insula and amygdala as well as interactions between the ACC and amygdala during the initial 1.2sec. Due to insufficient temporal resolution, the 2-sec model was unable to capture the rapid interactions that occurred in these critical emotion regions.
Of 66 interregional connections in the control group, 14 were between homologous regions, including amygdala, hippocampus, and anterior insula. Of 52 nonhomologous influences, 25 occurred in the first 1.2sec, with 9 from the right amygdala, 4 from the anterior insula, and 12 from the ACC. Of 26 interregional connections in the phobic group, the dorsal amygdala was responsible for 7 of the 8 influences between homologous regions (the 8th was again left to right anterior insula). Of 18 nonhomologous influences in the phobic group, 7 emanated from dorsal amygdala, 5 from anterior insula, and the remaining 6 from the ACC.
In the 400-ms model, there were 8 Bonferroni significant and 11 FDR significant path coefficients that were different between the groups. The majority of group differences (6 Bonferroni and 8 FDR) emanated from the hippocampi, and the remaining group differences were from subgenual ACC to ventral amygdala. For 7 of the 8 Bonferroni paths, phobics showed less inhibition than controls, as indicated by red in Figure 3C. Of particular interest is the lack of inhibition in phobics for the two paths from subgenual ACC to ventral amygdala, which had negative path coefficients in controls (−0.012 and −0.022 for 400ms and 2000ms, respectively) and were near zero for phobics (−0.0007 and −0.0005, respectively). The same pattern was observed for both hippocampus-amygdala paths (controls: −0.29 and −0.052 for 2400ms and 4000ms, respectively; phobics, 0.0006 and −0.0003, respectively).
To investigate if a particular context (i.e., trial type) drove the differences described above between snake phobic and nonphobic participants, we conducted context-dependent Granger causality analysis on the 400-ms data. Six trial types were individually investigated: disgust, fish, and snake, each with and without controllability. Of the 8 path coefficients that were different between the control group and the phobic group in the context-independent 400-ms model above, 4 were largely driven by a single trial type. The controllable snake trial type was responsible for both of the strongest path coefficients—left hippocampus to left ventral amygdala at 2400ms and right to left hippocampus at 3600ms (uncontrollable snake trial type for left hippocampus to left ventral amygdala at 4000ms; controllable disgust trial type for left hippocampus to right anterior insula at 2400ms; all ps<0.03).
For context-dependent Granger causality analysis implemented across all regions and interactions (Figs. S4–S6), the uncontrollable snake trial type had the most interregional paths showing Bonferroni-corrected group differences, all five of which emanated from the hippocampus, with phobics showing reduced influence for the four terminating in the amygdala.8
In addition to providing the direction of predictive influence between two regions (indicating which region does the influencing and which is being influenced), the implementation of Granger causality analysis (Chen et al., 2009) used here also provides the sign of the influence, which indicates whether the two regions change activity in the same (+) or opposite (−) direction. The 400-ms model showed a pattern of alternating positive and negative influences, as indicated by the color of the arrows in Figures 2 and and3.3. These positive and negative influences may indicate excitatory and inhibitory signaling pathways, respectively. However, other factors may contribute, such as an excitatory influence at one lag resulting in a negative influence at the next lag due to the increased activity being followed by a return to baseline.
This study showcases a new method for investigating effective connectivity of fMRI data at 400-ms resolution. Simulated data and real-world data presented here indicate the benefits of 400-ms resolution over native 2-sec resolution. Granger causality analysis of upsampled data from an anticipation paradigm using aversive video clips resulted in more influences between regions highlighted in prior research on anxiety and emotion than Granger causality analysis of data at native resolution. This improvement in the detection of effective connections was observed for healthy volunteers and snake phobics as well as for group differences. In addition, the method used here for upsampling to 400-ms resolution was superior to an alternative method in terms of preserving information about regional brain activity in the fMRI signal that is otherwise lost. Specifically, upsampling via interpolation at the voxel level using slice-timing correction to 400-ms increments resulted in a more dynamic time course than interpolation from the ROI average after conventional slice-timing correction. These data indicate that Granger causality analysis at native resolution overlooks many influences available in fMRI data and that increasing temporal resolution with interpolation at the voxel level is an effective method for investigating interregional influences at rates closer to those known to occur physiologically.
Granger causality at 400-ms resolution is an improvement over Granger causality at 2-sec resolution for multiple reasons. The simulated data at 400-ms resolution were closer to the ground truth than the data at 2-sec resolution, which was supported by the more dynamic time series for the real-world data at 400-ms resolution than at native resolution. In addition, the data at 400-ms resolution showed greater sensitivity in detecting effective connections than the data at 2-sec resolution. The increase in significant connections indicates an improvement, because the more rigorous Bonferroni correction for the 400-ms data ensures that the validity of the connections found are as valid as those for the 2-sec data. Finally, Granger causality at 400-ms resolution also provides greater temporal precision of identified connections. A path coefficient indicates the predictive influence of one region on another region at a given time lag. Multivariate regression solves for all of the different paths and time lags simultaneously, such that the values calculated for any path are corrected for the effects of all the other regions and time lags modeled. The 400-ms model has four time lags before reaching the first time point in the 2-sec model. As a result, influences can occur in the 400-ms model well before the first lag of the 2-sec model. This was evident in the current study, with Figure 3 illustrating numerous effective connections in the first four lags of the 400-ms model. It is likely that all interregional influences at even the first lag of the 2-sec model are the result of multiple synaptic connections that are hidden at this resolution.
A long-standing problem in fMRI research is the seeming limitations on temporal resolution resulting from the TR needed to acquire fMRI data for the whole brain. Highlighted in this study, voxel-level interpolation retains additional information about regional brain activity because multiple samples have been taken over the course of one TR before temporal blurring. Conversely, for ROI-level interpolation, the signals have been blurred and resampled once per TR before upsampling. Figure 1 illustrates the differences in the time courses for these two methods, indicating the increased information available for the upsampling method using voxel-level interpolation.
Another advantage to our upsampling method is that it allowed the investigation of context-dependent connectivity. The limited number of time points per condition and frequent breaks in the time series for the 2-sec data make context-dependent Granger causality analysis at 2-sec resolution unfeasible without extending the duration of each trial.1 However, upsampling to 400-ms resolution provided enough overlapping time points in each trial to allow for Granger causality analysis. Results here point to the potential of using context-dependent analysis for high temporal resolution Granger causality to elucidate the conditions that are driving connectivity patterns.
The upsampling method featured in the present study on effective connectivity is unlikely to confer advantages for many functional connectivity methods (e.g., psychophysiological interactions), because targets are individual voxels, which are only sampled once per TR and, therefore, cannot benefit from upsampling. To assess whether functional connectivity using ROIs rather than individual voxels as targets may benefit from upsampling, we conducted an exploratory analysis using a functional connectivity method examining correlations across the entire time series between each pair of 13 regions utilized in the current study. The results for the 400-ms model were essentially identical to the results for the 2-sec model, as would be expected given that one region's data points should show similar correlations with the same data points for another region regardless of whether those data points fall on the TR or on the 400-ms increments of the TR.
Although the focus here is on showcasing a new method for effective connectivity with fMRI data at 400-ms resolution, findings x also contribute to extant research on anxiety disorders. To our knowledge, there are no reports in the literature to date examining effective connectivity in specific phobia and only one for other anxiety disorders (Liao et al., 2010). Group differences were found only for effective connections emanating from the hippocampus and the subgenual ACC. Of particular interest, phobics showed less inhibition from the subgenual ACC to the amygdala than healthy controls, as would be expected based on previous findings for ACC-amygdala connectivity (Etkin et al., 2006; Johnstone et al., 2007; Phelps et al., 2004; Sarinopoulos et al., 2010; Urry et al., 2006). Interestingly, the phobic subjects also showed a similar reduction in inhibition from the hippocampus to the amygdala. Post hoc context-dependent analyses for these group differences revealed that the controllable snake trials were responsible for the two strongest group differences, both of which emanated from the hippocampus, suggesting an abnormality in phobics that is specific to the control of phobogenic stimuli. On the other hand, no single trial type was responsible for the group differences in subgenual ACC influences on the amygdala, suggesting a more general abnormality in this emotion regulation pathway (Etkin et al., 2010; Johnstone et al., 2007).
Granger causality analyses for each group separately revealed a wide array of effective connections among the amygdala, hippocampus, insula, and ACC, both intra- and interhemispherically. The extensive connectivity between the dorsal and ventral amygdalae and between the pregenual and supragenual ACC suggest that proximity may be a potential caveat, since both pairs of regions have shared surface boundaries. However, the pregenual and subgenual ACC are also in close proximity with shared surface boundaries, yet only one pathway was significant. Another important consideration is whether the presence of hemodynamic response variability across brain regions might complicate network detection (Chang et al., 2008; David et al., 2008; Deshpande et al., 2010b; Friston, 2009a), particularly for the within-group analyses. For between-group analyses, any hemodynamic response variability across brain regions would be removed when comparing the two groups, unless there were also reliable group differences in hemodynamic response variability. Future research in the area could apply arterial spin labeling (Detre et al., 1992; Wong et al., 1997) and breath holding (Birn et al., 2008; Chang et al., 2008), which are both promising approaches for determining subject-specific regional hemodynamic response delays that might inform contributing factors to significant interregional influences observed within and between groups. Additionally, between-group analyses have the advantage of being entirely resistant to artifact from deterministic data transformation, such as upsampling, which is a potential confound for within-group analyses.
A potential confound for applying effective connectivity methods to upsampled data is that auto-correlation between time points is inflated. For Granger causality, auto-correlation is modeled with auto-loops and, therefore, is separated from interregional connections (Table 1), although often ignored in other reports using Granger causality. Moreover, the simulated data further attest to the validity of the interregional connections found, in that the strength of the Granger causality weights increased as the ground truth connectivity of the simulated data increased. Consequently, the larger number of significant interregional pathways found here for 400-ms resolution than for 2-sec resolution cannot be attributed to differences in auto-correlation. In addition, it is highly unlikely that significant connections here are false positives because the stringent Bonferroni correction procedure utilized is the most conservative method for preventing false positives.
One limitation of the current study is that neural events occur at time scales faster than 400ms. Future investigations could directly test whether models at finer resolutions (that may detect phase shifts less than 400ms between regions) are superior at detecting regional interactions, even with the greater stringency of statistical correction procedures for comparing to other models encompassing an equivalent time frame (see “Data acquisition, processing, and upsampling”). The present report is a promising start to this line of work, with the 400-ms model indicating improved detection compared to the standard 2-sec model for time frames of 4 and 6sec, respectively, despite a statistical threshold that was 333% more stringent. Moreover, the upsampling approach here is also relevant for the recently developed structural vector autoregression method of effective connectivity for further addressing temporal resolution in fMRI data (Chen et al., 2011).
Granger causality analysis of fMRI data that was upsampled to 400-ms resolution using a new method identified many causal network connections that were overlooked by analysis at native resolution. The upsampling method introduced here capitalizes on interleaved slice acquisition and slice-timing correction to identify rapid causal influences that occurred in less time than the TR duration. The benefits of this method for improving temporal resolution in fMRI were observed for simulated data and real-world data, and for both context-independent and context-dependent Granger causality. Effective connections were detected at far higher levels for 400-ms resolution than 2-sec (native) resolution in healthy volunteers and a clinical population with snake phobia, as well as for group differences comparing the phobics to healthy subjects. Taken together, data here indicate the promise of using fMRI to investigate causal influences between brain regions at rates closer to those known to occur physiologically.
1For the case of 30-slice acquisition, 29 slices are interpolated to 1 slice for each time increment. For data at native resolution, the time of only 1 slice is used for interpolating the other 29 slices, whereas for data at upsampled resolution using our method of correcting to multiple slices, the time of each of those slices is used for interpolating the other 29 slices, again resulting in the proportion of 29/30 or 97%.
2Informed consent was obtained from all participants before the experiment, in accordance with rules set by the University of Wisconsin School of Medicine and Public Health Institutional Review Board.
3An important consideration when upsampling a signal is the possibility of aliasing frequencies above the Nyquist frequency. Since each voxel in the current study was sampled and interpolated from 0.5-Hz data (corresponding to TR=2), only frequencies up to the Nyquist frequency of 0.25Hz can be accurately represented through Fourier interpolation. Frequencies higher than the Nyquist are not included in Equation S1 or S2 of Supplementary Methods and could be mis-assigned to one of the m frequencies used in those equations, a phenomenon known as aliasing (Fig. S3A). However, since the slices are acquired interleaved, many frequencies higher than the Nyquist will be represented by opposite signs in adjacent slices and, thus, average to zero, canceling out those high frequencies (Fig. S3B). Therefore, the interleaved acquisition helps to prevent frequencies beyond the Nyquist from being included in the signal. Of note, frequencies more than approximately double the Nyquist could still alias with both slices having the same sign. However, frequencies with such rapid oscillations are not observed or expected in fMRI signal, since the hemodynamic response (HDR) is a relatively slow, low-frequency signal. Also of relevance for effective connectivity, a causal relationship between two brain regions can only be inferred if they have multiple frequencies that show similar differences in phase, such that a single aliased frequency should have little effect on the model.
41dGC.R uses four criteria (AIC, Final Prediction Error, Hannan-Quinn, Schwarz-Bayes Criterion) to guide how many lags should be modeled. The AIC recommendation was 3 lags for the 2-sec model and 10 lags for the 400-ms model.
5There are three advantages of 1dGC.R not shared by all Granger causality methods: (1) the sign of the influence is maintained, (2) the influence from one region to another is tested for each lag independently, (3) the two-stage analysis noted in the text that performs Granger causality on individual subject data followed by second-level analysis to determine statistically consistent influences across subjects. As such, the Granger causality method used here produces statistics that reflect consistency across subjects in direction and magnitude for each time lag.
6The number of data points needed for a model to converge is dependent on the data (e.g., signal-to-noise ratio, effect size). There are 132 total trials (33 trials in each of four runs) averaging 22sec per trial, with 2-sec TRs. Hence, for a simple correlation analysis with no lags, there are 1452 time points (132*22/2) for each subject. For Granger causality, the number of time points is reduced by the number of lags multiplied by the number of breaks in the time series (including the end). The 2-sec model with three lags and four runs (breaks) has 1440 data points included in the analysis (1452–3*4). To analyze context-dependent Granger causality, the appropriate trial types must be separated. As a result, the number of breaks increases from 4 to 132 (corresponding to the number of runs and trials, respectively), reducing the number of data points to 1056 (1452–3*132). In addition, these 1056 data points are distributed over six different conditions, so the model for each condition has only 176 data points (as point of comparison, each condition for the 400-ms model has 885 data points). It is possible that a model may converge with 176 data points for some subjects under some conditions. However, to do a group-level analysis, the model must converge for all subjects and conditions, which happened for the 400-sec model, but not for the 2-sec model.
7Two primary multiple comparison correction methods—Bonferroni and false discovery rate (FDR)—were presented in the results, because they address different objectives. Bonferroni at p<0.05 is a highly conservative method that aims for 5% chance that any of the significant paths are false positives. FDR at p<0.05 is less conservative in that it aims for 5% of the significant paths to be false positives.
8Of note, differences between groups were not exclusively driven by a lack of connectivity in the phobic group. For example, the FDR-corrected group difference for the uncontrollable snake trial type from right ventral amygdala to left anterior insula at 2000ms, an influence was observed in phobic subjects (0.20), but not controls (−0.0041).
This work was supported by the National Institute of Mental Health [R01-MH74847, K08-MH63984, and K02-MH082130 to JBN]; and by a core grant to the Waisman Center from the National Institute of Child Health and Human Development [P30-HD03352].
We gratefully acknowledge Andrew Alexander, Michael Anderle, Danielle Green, Ron Fisher, Donald McLaren, Desmond Oathes, John Ollinger, Adrian Pederson, Hillary Schaefer, Allison Schaus, and Nathan Vack for their contributions to this project.
No competing financial interests exist.