|Home | About | Journals | Submit | Contact Us | Français|
Correlation and causality metrics can be applied to blood-oxygen level dependent (BOLD) signal time series in order to infer neural synchrony and directions of information flow from fMRI data. However, the BOLD signal reflects both the underlying neural activity and the vascular response, the latter of which is governed by local vasomotor physiology. The presence of potential vascular latency differences thus poses a confound in the detection of neural synchrony as well as inferences about the causality of neural processes. In the present study, we investigate the use of a breath holding (BH) task for characterizing and correcting for voxel-wise neurovascular latency differences across the whole brain. We demonstrate that BH yields reliable measurements of relative timing differences between voxels, and further show that a BH-derived correction can impact both functional connectivity maps of the resting-state default-mode network and activation maps of an event-related working memory (WM) task.
There is a growing interest in using fMRI to study distributed networks of interacting regions across the whole brain (Marrelec et al., 2006; Rogers et al., 2007; Varela et al., 2001). Currently, time series modeling techniques are applied to regional BOLD signals in order to infer synchrony and directional influence among neural populations, both within specific tasks (see, e.g., (Marrelec et al., 2006; Rogers et al., 2007) for review) and, more recently, during task-free resting state (see (Fox and Raichle, 2007) for review).
A fundamental assumption of most fMRI correlation and causality studies is that the relative timing of BOLD signals from different regions of cortex reflects the timing of their underlying neural activity. However, the BOLD signal results from a coupling between neuronal activity and vascular responses (Bandettini et al., 1992; Kwong et al., 1992; Logothetis et al., 2001; Ogawa et al., 1992). Thus, differences in the timing (latency) of vascular responses between brain regions will affect the relative timing of the BOLD signals, thereby confounding the ability to infer neural connectivity and causality. The dynamics, location, and magnitude of the BOLD signal are heavily influenced by the vasculature in each voxel; large vessel effects can cause delays of up to 4 sec relative to capillary effects (Bandettini, 1999; Lee et al., 1995). Widely disparate hemodynamic response functions (HRFs) have been found across the cortex(Bandettini, 1999; Buckner et al., 1998; Handwerker et al., 2004; Robson et al., 1998; Schacter et al., 1997; Thomason et al., 2005), and even adjacent voxels can exhibit latency differences of up to 2 sec, which are most likely of vascular origin (Miezin et al., 2000). Latency differences may become more pronounced at higher spatial resolutions due to greater variability of partial volume effects from vascular compartments. Therefore, it is critical to model and correct for non-neural latency differences in the BOLD signal before attempting to quantify inter-regional correlations or causality.
While several researchers have quantified HRFs and timing delays in focal (sensory) regions prior to modeling their interactions (e.g. (Chen and Desmond, 2005; Menon et al., 1998; Miezin et al., 2000)), the vascular latency confound is largely overlooked in studies of whole-brain functional connectivity. One reason is due to the difficulty of characterizing HRFs across the whole brain, as most tasks used to probe HRFs activate only a small set of regions robustly and with adequate signal-to-noise ratio, and virtually all such HRF measurements have been limited to sensory regions. Even if such measurements could be made with sufficient reliability and accuracy, it is not clear that they generalize to other regions because of known regional differences in the vascular system. A second reason is that whole-brain, resting state studies are less amenable to analyses such as condition-dependent modulation of connectivity/causality (e.g. (Roebroeck et al., 2005; Sun et al., 2007)) which essentially subtracts away the common vasomotor component of the HRF delay and thus circumvents the need for latency correction.
In the present study, we investigated whether a hypercapnic challenge (breath holding, or BH) can be used to quantify vascular latency differences across the whole brain. BH causes local up-regulation of blood flow when reduced perfusion, in the presence of ongoing baseline metabolism, leads to vasodilation of capillaries that is characteristic of local reactivity (Corfield et al., 2001; Kastrup et al., 1999a; Kastrup et al., 1999b; Kastrup et al., 1999c; Kastrup et al., 1998; Liu et al., 2002; Nakada et al., 2001). Since BH modulates cerebral blood flow to all vascularized brain regions without an accompanying change in CMRO2 (Kastrup et al., 1999a), it is a simple and robust method for assessing regional vascular reactivity properties, uncoupled from neural activation (Bandettini and Wong, 1997; Cohen et al., 2004; Thomason et al., 2005; Thomason et al., 2007). We hypothesized that the latency of a voxel in the BH task would reflect the vascular component of its HRF latency, and that adjusting its time series accordingly would reduce confounds in further analysis.
Although previous studies have found an empirical relationship between BH and task-activation magnitudes (Handwerker et al., 2007; Thomason et al., 2007), it is not known whether BH timing delays will be directly proportional to activation-induced vascular delays. The BH latency may incorporate the transit time of blood flow from major artieries in addition to the intrinsic reactivity of local capillary beds, and it is possible that the former may dominate the latency measurement in BH. In an analysis of time-to-peak and coherence phase, (Handwerker et al., 2007) found that latency in a BH task was not reliably predictive of latency in a visuomotor saccade task, though occasional significant correlations were observed; however, the relationship may have been weakened due to behavioral variability in motor response time during the saccade task. In the current study, we investigated the agreement between latency measured by BH and that from a paradigm designed to evoke simultaneous activation of primary sensory (auditory, visual, and motor) regions, with the expectation that a correspondence between latency values across the 2 tasks would provide evidence that BH latency primarily reflects local vasomotor reactivity.
As described previously, regional differences in latency would impact functional connectivity maps of networks such as the default mode network (DMN). The DMN encompasses a set of regions that exhibit low-frequency correlated signals in task-free resting state (Greicius and Menon, 2004; Raichle et al., 2001), and collectively down-regulate during a wide range of cognitively demanding tasks (Binder et al., 1999; Gusnard et al., 2001; Mazoyer et al., 2001; McKiernan et al., 2003; Shulman et al., 1997). DMN connectivity has been shown to vary across groups (Baliki et al., 2008; Damoiseaux et al., 2007; Garrity et al., 2007; Greicius et al., 2007; Greicius et al., 2004; Uddin et al., 2008) and cognitive states (Esposito et al., 2006; Waites et al., 2005), and an increasing number of studies seek to make inferences about behavior (Clare Kelly et al., 2008; Daselaar et al., 2004; Hampson et al., 2006) and dysfunction (Garrity et al., 2007; Greicius et al., 2007; Greicius et al., 2004; Uddin et al., 2008) from measurements of DMN connectivity. Removal of non-neural latency differences is important for obtaining interpretable, quantitative results, particularly for between-group studies in which hemodynamics may pose a major confound.
Latency is also an important consideration in detecting task activation, particularly when event-related designs are employed. Analysis of event-related designs is sensitive to time shifts between the actual and modeled voxel responses, and methods to estimate and compensate for potential delays within the framework of the general linear model have been proposed (Calhoun et al., 2004; Friston et al., 1998; Liao et al., 2002; Worsley and Taylor, 2006). In the present study, we also examined the effect of a BH latency correction on activation maps in an event-related working memory (WM) task.
Lastly, to illustrate the potential impact of BH latency correction on causality analysis, we examined correction-induced changes in Granger causality among a set of regions in the WM task for one subject.
Participants included 10 healthy adults (4 female) between the ages of 22 and 42 (mean age = 27.6), including 1 left-handed male. All subjects provided written, informed consent, and all protocols were approved by the Stanford Institutional Review Board.
In the BH scan, subjects were monitored by a respiratory belt (TSD 201, Biopac Systems, Santa Barbara, CA) placed snugly around their upper thorax. The belt’s electrical conductance is nominally proportional to the belt circumference and thus varied as the subject breathed. The belt was energized by a 5-V power supply through a 24-Ω resistor and connected to a custom-made analog to digital (A/D) converter, which in turn was connected to the stimulus presentation computer as a parallel device. During a preliminary setup step, subjects were instructed to first take in a deep breath, during which the value was recorded, and then let their breath out while a second recording was made. The subject’s range of breathing during this calibration was used to obtain a target inspiration depth for breath holding (‘target zone’, which was 90% of the maximum inspiration), so that the inspiration depth was accurately monitored during the BH task described below.
During all functional scans, subjects’ heart rate and respiration were monitored at 40 samples/sec using the scanner’s built in photoplethysmograph placed on the left hand index finger and a pneumatic respiratory belt strapped around the upper abdomen, respectively.
Subjects underwent 4 functional scans: (1) a BH task, (2) a sensory-motor (SM) task, (3) a working memory (WM) task, and (4) resting state (no task). All stimuli were programmed and presented in E-prime (http://www.pstnet.com).
Subjects performed 16 repetitions of a block that included a “normal breathing” phase and a “breath hold” phase. A schematic of one block is shown in Figure 1. The timing and depth of breathing were cued by auditory tones. A low-pitched tone cued the subject to take a normal breath; after 3 normal breaths, a medium-pitched tone (‘get ready’) was presented, upon which the subject would take a normal breath but anticipate the following cue, which was a ‘deep breath’ cue. A high-pitched tone was subsequently presented, cueing the subject to take a deep breath (reaching their target zone within 2 sec; a rapid double-beep sounded when the target zone was attained) and hold it for 11 sec. A second high tone sounded to signify the end of the 11 sec period. The “normal breathing” tones were separated by 3.5 sec, as pilot testing had indicated that this was a comfortable breathing rate for the individuals tested. After the final (16th) block, an additional 28 sec block of “normal breathing” was appended in order to capture the full hemodynamic response to the final “breath hold” phase.
The auditory cues undoubtedly introduce a small amount of activation in the auditory cortex; however, task-correlated activation was minimal (as determined by pilot scanning in which subjects only listened to the stimuli without following the breathing cues). Furthermore, in a series of pilot runs, we observed a strong correspondence between latency measurements in the sensory areas derived using the auditory stimulus and those derived using a visual stimulus, and when varying the breathing rate and duration of the breath hold phase.
The SM task consisted of alternating “on” and “off” blocks of varying length, with block lengths and onsets chosen to increase efficiency and detection power beyond those of conventional periodic block and event-related designs. “On” blocks consisted of simultaneous visual, auditory, and sensorimotor stimulation. Visual stimuli consisted of a circular checkerboard of alternating black and white contrast, reversing at 3 Hz. Auditory stimuli consisted of a randomized binaural tone sequence (3 Hz, synchronized with the visual stimuli; tone duration = 166 ms, interval between tones = 167 ms). Passive motor stimulation was delivered by pneumatically-driven plungers arranged in left and right hand gloves so as to contact and push all 10 fingers up and down in a randomized pattern at 3 Hz, in synchrony with the checkerboard reversals and auditory tones. The total scan time was 5.5 min. The use of such a passive task minimized potential behavioral differences across subjects in performing finger apposition.
Subjects performed an event-related Sternberg verbal WM task (Sternberg, 1969), consisting of 3 back-to-back 5 sec trials (0.5 sec encoding, 3 sec maintenance, 1.5 s probe) presented with a mean ISI of 60 sec, with a total scan duration of 10 min. Encoding stimuli consisted of 4 uppercase letters in a cross-like configuration around the center of the screen, and the probe stimulus was a single lowercase letter presented at the center of the screen. In addition to examining the effect of the BH latency correction on activation maps, this task was used to localize subject-specific seed regions for deriving the DMN in the resting state scan (described below). Data for the WM task were unavailable for 2 subjects; one subject was scanned before the WM task design was implemented, and another subject’s imaging data file was corrupted.
In the 8-min resting state scan, subjects were instructed to relax and close their eyes while remaining awake.
Magnetic resonance imaging was performed on a 3.0-T whole-body scanner (Signa, rev 12M5, GE Healthcare Systems, Milwaukee, WI) using a custom quadrature birdcage head coil. Head movement was minimized with a bite bar. Thirty oblique axial slices were obtained parallel to the AC-PC with 4-mm slice thickness, 1-mm skip. High-resolution T2-weighted fast spin echo structural images (TR = 3000 ms, TE = 68 ms, ETL = 12, FOV = 22 cm, matrix 192 × 256) were acquired for anatomical reference. A T2*-sensitive gradient echo spiral in/out pulse sequence (Glover and Lai, 1998; Glover and Law, 2001) was used for functional imaging (TR = 2000 ms, TE = 30 ms, flip angle = 77°, matrix 64 × 64, same slice prescription as the anatomic images). A high-order shimming procedure was used to reduce Bo heterogeneity prior to the functional scans (Kim et al., 2002). Importantly, a frequency navigation correction was employed during reconstruction of each image to eliminate blurring from breathing-induced changes in magnetic field; no bulk misregistration occurs in spiral imaging from off-resonance (Pfeuffer et al., 2002).
No motion coregistration was performed, as (1) motion was expected to be minimal because a bite-bar was used, (2) coregistration causes unintended smoothing across voxels, which would interfere with our voxel-wise analysis, and (3) estimation of coregistration parameters can be biased by activation in tasks (such as BH) that evoke large, widespread signal changes (Freire and Mangin, 2001). Therefore, it was important to verify that intra- and inter-scan motion was minimal. We calculated the maximum peak-to-peak excursion, RMS fluctuation, and task-correlated motion for the 3 translational and 3 rotational motion parameter time series within each scan (Thomason and Glover, 2008). Because the latency values in the BH scan were applied in a voxel-wise fashion to the SM, WM, and rest scans, we also report the total drift across the BH scan and each of the 3 other scans. Summary statistics are reported as the maximum of these values over the 6 axes of motion.
Functional images were pre-processed using custom C and Matlab routines. Pre-processing included slice-timing correction, spatial smoothing with a 3D Gaussian kernel (FWHM=5mm), and removal of linear and quadratic temporal trends. Spatial normalization to a standard template was not performed as a pre-processing step; all computations were done in the original subject-space, and voxels were maintained in their original dimensions (3.4375mm × 3.4375mm × 4mm). The first 7 temporal frames were discarded to allow the MR signal to equilibrate. In the BH task, temporal frames corresponding to the first 2 blocks of the design (32 frames) were discarded because subjects performed less consistently at the beginning of the task. Resting state data were further corrected for cardiac and respiratory artifacts using RETROICOR (Glover et al., 2000) and correction for respiratory variation, which has been shown to account for a substantial amount of the variance of resting state global fluctuations (Birn et al., 2006) Briefly, the correction involved sampling the standard deviation of the respiratory waveform over every TR-interval (2 sec), convolving it with a respiratory variation impulse response function (RRF(t) = 0.6t2.1et/1.6 − 0.0023t3.5e−t/4.25) (Birn et al., 2007) and regressing it out of each voxel time series using cross-correlaton at shifts from −8 sec to 8 sec, in intervals of 1 sec.
Voxel-wise hemodynamic latencies were measured using linear cross-correlation between each voxel time series x(t) and a reference time series y(t). The latency was defined as the time-lag yielding the maximum of the cross-correlation function:
Both x(t) and y(t) were first upsampled to a resolution of 100 msec using Fourier interpolation, and τ was permitted to range between ±4 sec.
For the SM task, the reference time series y(t) consisted of the binary design waveform convolved with a canonical HRF (Glover, 1999). For the BH task, y(t) consisted of the average time series across voxels in the whole brain which exceeded a minimum correlation of r>0.25 with the canonical BH task regressor (which typically included most gray matter voxels). Using a brain-derived reference in the BH task, rather than an external reference, accounts for the variations in a subject’s performance in the BH task that occur despite the auditory cues and feedback. It is also appropriate for our goal of quantifying relative latency differences across voxels, rather than absolute latency. Relative latency differences were found to be preserved for various alternative functions that were tried (e.g. using an average time series from sensory regions activated in the SM task).
To assess the validity of BH as a method for measuring hemodynamic latency, we compared latency measurements across the BH and SM tasks by computing the Pearson correlation coefficient between voxel-wise latencies in the 2 tasks. We expected that strong positive correlations between latency in the BH and SM tasks would indicate that BH yields meaningful measurements of local vasoreactivity, and, by extension, that BH latency measurements would also apply to brain regions outside of these sensory cortices.
The precision of the cross-correlation method in estimating latency depends on the contrast-to-noise ratio (CNR) of the voxel time series. To quantify the latency precision of our BH and SM tasks, simulated noise was added to the HRF-convolved stimulus waveform of each task at CNRs varying from 0.1–2, and the standard deviation of the lag estimate was computed over 1000 independent iterations. We used both Gaussian white noise and AR(1) noise, as fMRI time series have been shown to exhibit temporal autocorrelations. An AR(1) parameter of 0.5 was used, as it reflected the mean AR parameter estimate of time series from a number of gray-matter voxels sampled from resting state data. CNR was transformed to Pearson correlation coefficients using the relationship
To correct for neurovascular latency differences in the resting state, the time series of a voxel in the resting state scan was first upsampled to a resolution of 100 ms using Fourier interpolation, and then time-shifted by the corresponding latency value obtained from the BH task. The correction was applied to all voxels which exceeded r>0.6 (t(222)=11) in the BH task, and was performed on a single-subject basis.
Simply shifting the time series by the latency value measured in the BH task is sufficient for analyses in which only the relative latency differences among voxels matters, such as in functional connectivity and causality. However, when assessing task activation using a GLM, it is critical to have an absolute reference point – i.e., a region with known hemodynamic latency relative to stimulus onset – from which to perform all latency shifts; such an absolute reference is lacking from the BH latency maps due to the use of a brain-derived waveform for the cross-correlation analysis. One can define a reference point on a subject-by-subject basis by measuring the HRF latency in an ROI using a task with robust activation (e.g. a sensory task) and offsetting the BH latency map accordingly. To correct for latency differences in the WM task, for each subject we used the intercept of the least-squares fit between voxel-wise latencies in the BH and SM task to adjust the offset of the BH latencies. After upsampling to a resolution of 100 ms, the time series of suprathreshold (r>0.6 in the BH task) voxels in the WM task data were then shifted by their corresponding adjusted BH latency values.
Functional connectivity maps of the DMN were compared before and after the latency correction. At the individual subject level, a DMN functional connectivity map was obtained by extracting the average time series from a seed ROI in the posterior cingulate cortex (PCC)/precuneus, a central node of the DMN (Greicius et al., 2003; Greicius and Menon, 2004; Raichle et al., 2001). The seed ROI was constructed by forming an 8 mm sphere centered at the peak deactivated (rest > experimental condition) coordinate in the WM task such that the entire sphere fell within a template anatomical PCC/precuneus (MarsBar AAL atlas; http://www.sourceforge.org/marsbar). For the 2 subjects in which data from the WM scan was unavailable, the coordinate of peak deactivation in the SM task (which also elicited PCC deactivation) was used instead. The final seed ROI was comprised of voxels in the sphere which exceeded the BH correction threshold (r>0.6), so that the same seed ROI could be used for both the uncorrected and latency-corrected analyses. We verified that at least 80% of the original 8mm sphere was retained in this step, so that ROI sizes were comparable across subjects. The time series from the seed ROI was then correlated (at zero lag) against that from every voxel in the brain, generating a map of correlation coefficients.
To test for significant connectivity changes in an individual subject, a difference map between the Fisher-Z transformed corrected and uncorrected maps was normalized to a Z statistic map ( , here with n = 233 frames (Dowdy and Wearden, 1991)). To examine whether any common changes in network connectivity existed across the group, single-subject difference maps (between Fisher-Z transformed corrected and uncorrected connectivity maps) were normalized to a standard EPI template and entered into a random-effects group-level analysis using SPM5 (http://www.fil.ion.ucl.ac.uk/spm).
WM activation maps before and after the latency correction were obtained by correlating the time series of each voxel with a canonical task waveform. Correlation coefficients were transformed into Fisher Z values, and the single-subject difference maps between the corrected and uncorrected maps were normalized to a standard EPI template and entered into a random-effects group-level analysis using SPM5.
Granger causality analysis was performed in the WM task for one subject. A subset of regions known to show increased activation or deactivation in WM were selected for the analysis; activated regions included dorso-lateral prefrontal cortex (DLPFC), insula, posterior/superior parietal cortex, and anterior cingulate cortex (ACC), and deactivated regions included the PCC and ventro-medial prefrontal cortex (VMPFC) (Smith et al., 1998). ROIs for the analysis were defined by taking the intersection of the above anatomical ROIs with the subject-specific task activation map (liberally thresholded at |r|>0.2), and selecting the voxels which exceeded the threshold for latency correction (r>0.6 in the BH task). Granger causality was implemented with the aid of the Causal Connectivity Toolbox (Seth, 2005), and significant (p<0.05, Bonferroni corrected) multivariate causal connections were compared before and after latency correction.
Subjects had minimal overall and task-correlated motion in the BH, SM, and WM tasks (Table 1). Over the combined scans (BH+SM, BH+WM, and BH+Rest), a drift of less than one voxel was observed for all but one subject, who showed a 6.9 mm drift across the BH and WM tasks combined. It was found that results from this subject were consistent with the others, so the data were not excluded. A second subject had excessive motion only during the SM task (he found the stimulus tones amusing and laughed throughout the task) and so was excluded from all analyses.
Figure 2 shows the standard deviation of the latency estimate versus correlation coefficient for the BH and SM tasks, using simulated white noise and AR(1) noise. To choose a suitable correlation threshold, we examined typical activation maps. The BH task elicited robust global activation at r>0.6 (p<1e-18, Bonferroni correction; Figure 3a), as assessed using cross-correlation with a [negative] canonical HRF convolved with the task design. The SM task induced strong (r>0.6, p<1e-12, Bonferroni correction) responses in the visual, auditory, and sensorimotor cortices (Figure 3b). With reference to Figure 2, the distributions of cross-correlation values in both tasks suggest that latency can be estimated to a worst-case variability of ~500 ms in activated voxels, using liberal thresholds such as r>0.6.
In a pilot study of the BH task design, latency measurements were replicable across the preferred (auditory cueing) BH task and a variant of the BH task that used visual cues. We obtained a correlation of r=0.77 (p<1e-62) between the voxel-wise latency measurements across the 2 tasks.
Figure 4 shows latency maps from the BH task and SM tasks, along with the distribution of latency values in both tasks, for 2 different subjects. Latency tended to exhibit smooth variations across space, but the degree of overall heterogeneity varied across subjects. While it is reasonable to expect that the spatial pattern of hemodynamic latency might vary considerably across individuals, we noticed that several features of the latency maps (notably, delayed responses in the cerebellum as well as in major vessels, such as the sagittal sinus) were common across the group. A group-average latency map, constructed from normalizing the subject-specific latency maps (each thresholded at a cross-correlation of r>0.6; t(222)=11) to a standard EPI template, is shown in Figure 5.
We obtained significant positive correlations between voxel-wise latency measurements in the (within-subject) BH and SM tasks (Table 2). Voxels included in the comparison were the intersection of those surviving a threshold of r>0.8 in the BH task and r>0.6 in the SM task. A higher threshold was selected for the BH task since the distribution of cross-correlation values was skewed toward higher values compared to those in the SM task; if the SM were thresholded as stringently, there would be too few suprathreshold voxels to conduct an adequate comparison. In addition, voxels falling at the boundaries of the interval used to assess latency (−4sec, 4 sec) were excluded, as any voxels with latency outside this interval in both tasks would be rounded to (−)4 sec and thus inflate the correlation. To investigate whether the reported BH v. SM task correlations were sensitive to the specific choice of thresholds, we repeated the analysis sweeping over a range of T-value thresholds (with T chosen identically for BH and SM), and did not find significant variation in correlation for thresholds yielding approximately 100 to 2000 voxels.
For the majority of subjects, the mode of the latency distribution (usually at or near zero lag) was often much higher than the remainder of the distribution (see, e.g., the histograms in Figure 4). If a large, common set of voxels were exhibiting zero relative lag in both the SM and BH tasks, it might bias the significance of the between-task correlation. Therefore, we repeated the BH v. SM task correlations after excluding the voxels in the mode of each task. As shown in Table 2A, correlations remained highly significant.
It was possible that the observed correlations were driven primarily by global properties of the latency maps (e.g. overall regional latency differences in the visual, compared to auditory, cortices). To examine whether the latency correspondence between BH and SM also held on a more local basis, we correlated voxel-wise latency values in the visual, auditory, and motor cortices separately. We observed a highly significant correspondence in all 3 regions in most subjects (Table 2B).
A perfect correspondence between BH and activation-induced responses would be reflected not only in strong positive correlations, but also in obtaining a regression-line slope of 1. For the sets of voxels analyzed at the reported thresholds, the majority of slopes (in plots of SM (y-axis) against BH (x-axis)) was less than 1 (Table 2). Yet, further analysis showed that (1) across the 26 sets of voxels (pooled over 9 subjects) shown in Table 2B, there was a significant relationship (p<1e-6) between the strength of the latency correlation and the proximity of the slope to 1 (Figure 6A), and (2) we observed that the distribution of slopes became closer to 1 as higher thresholds were applied to select voxels in the BH task (Figure 6B). Together, these provide further evidence for BH latency as a valid model, and suggest that the model becomes more accurate as the precision of the BH measurement increases, in accordance with our simulations (Figure 2).
Functional connectivity maps of the DMN before and after latency correction are shown for 2 subjects in Figure 7. At the individual subject level, correction resulted in significant connectivity changes, though with a modest spatial extent (Table 3); the percentage of voxels exhibiting significant connectivity changes ranged from 0–15% across subjects (mean=3.6%, SE=0.01%, N=9) at the p<0.05 uncorrected level, and from 0–7% (mean=0.01%, SE=0.007%, N=9) at the FDR=0.05 level. At the group level, significantly higher post-correction connectivity was confined to a small region of the cerebro-spinal fluid, which is most likely artifactual. More importantly, significantly lower (p<0.05) post-correction connectivity was widespread and localized primarily within gray matter (Figure 8A).
Significant group-level changes in the WM activation maps due to the latency correction are shown in Figure 9. Both increaes and decreases in activation were localized to areas known to be involved in WM (Cabeza and Nyberg, 2000; Smith and Jonides, 1999). Increased activation was observed bilaterally in the occipital cortex, caudate/putamen, and postcentral gyri, while decreased activation was observed bilaterally in the mid/inferior temporal cortex, angular gyri, ventro- and dorso-lateral prefrontal cortex, and cerebellum.
Figure 10a illustrates changes in the multivariate Granger causality analysis due to latency correction, for one subject. Before latency correction, significant causality was observed from the right insula to the ACC, left DLPFC, and left parietal cortex; after the correction, these 3 connections were preserved and several others became significant, such as bi-directional connections between the left and right parietal cortices as well as connections among the parietal and prefrontal areas and from the PCC.
We also examined the results of a pairwise Granger causality analysis on the BH data of this subject, using the same ROIs (Figure 10b). We choose to show pairwise, rather than multivariate, results for the BH because a pairwise analysis is more revealing of the directionality among all regions (multivariate analysis will omit a connection between 2 regions when the variance is better explained by common influence from a third). One observation is that the BH causality map shows a causal connection from right insula to the DLPFC, parietal, VMPFC, and ACC, suggesting that the finding of causality between these regions in the WM task from the may be explained by vasoreactivity; however, the persistence of these arcs from the right insula after the latency correction increases the likelihood that it may also be neurally-driven.
We have proposed the use of a hypercapnic challenge (BH task) as a method of measuring and correcting for relative hemodynamic latency differences across the whole brain. Such a task allows the interrogation of most gray-matter regions of the brain (Figure 3) and is thus not limited to specific regions, as is the more typically employed visual (or other sensory) task. Our results indicate that BH is a robust method for assessing non-neural, vasoreactivity-based latency differences across individual voxels. While spatial variability of the HRF shape beyond latency differences may still influence fMRI analyses, we propose that correcting for vasoreactivity latency confounds in an unbiased manner using a BH task will increase the validity of whole-brain functional connectivity and causality analyses.
The present study assumes that latency measures derived from a BH task are uncoupled from neural activation associated with either the performance of the BH task or the mild hypercapnia itself. Studies with concurrent BOLD and CBF measurements using the Davis model (Davis et al., 1998) have shown that CMRO2 changes during a BH task are small (Kastrup et al., 1999a), thus indicating that neural processes are unlikely to contribute greatly. BOLD signal magnitudes measured during BH were demonstrated to correlate well with those of hypercapnic challenges induced by CO2 inhalation, which requires no active participation (Kastrup et al., 2001); furthermore, Thomason et al. found that different forms of subject feedback used to control depth of inspiration during a BH task did not alter the BH response (Thomason and Glover, 2008), again suggesting that intentional BH control mechanisms did not significantly contribute to the overall BOLD signals. Finally, cognitive processes generally induce much smaller BOLD signals than are observed in BH (which are on the order of 2–4% or larger, even in motor and thalamic regions (Kastrup et al., 1999b; Thomason et al., 2005)) and thus would not add significantly to the vasoreactive response.
Previous studies have addressed latency confounds in whole-brain analysis by using cross-correlation rather than zero-lag correlation to build connectivity and activation maps (Calhoun et al., 2003). Cross-correlation effectively reduces the sensitivity of detection power to minor latency shifts, for voxels with adequate CNR. However, cross-correlation fails to distinguish between neural and hemodynamic delays, and furthermore uses the same dataset to both estimate and correct for latency; hence, the resulting maps are biased in an unknown way. On the other hand, BH elicits a strong BOLD signal change in all vascularized regions (e.g. Figure 3), and is furthermore unbiased since the correction is applied prior to functional analysis and the correction parameters are acquired from a separate scan. BH is also simple and practical to implement; it is non-invasive, comfortable for subjects to perform, and involves a relatively short additional scan (<8 min).
Performing a latency correction using BH relies on the assumption that the timing of local vasomotor responses to a BH task is proportional to those of neural activation-induced responses. We examined the validity of this assumption by correlating voxel-wise BH latency measurements with those of a task designed to elicit simultaneous activation of sensory regions. Our results showed strong positive correlations, with best-fit slopes approaching unity in the limit of more precise BH measurements, thus supporting the hypothesis that BH latency is reflective of local vasoreactivity delays and that the transit time of blood in the major arteries may be ignored.
Because regions of the cerebellum tended to exhibit high latency, we considered the possibility that the cerebellum may be more susceptible to transit time effects than other regions. However, in a subject with high cerebellar latency in the BH task, we also observed large time lags in a cross-correlation between the PCC and the cerebellum in the (uncorrected) resting state data (Figure 11). This finding helps suggest that the origin of the cerebellar delay is either vasoreactivity or another factor, such as respiratory effects.
Latency correction yielded differences in subject-specific DMN maps. Most differences were minor, perhaps owing to the fact that connectivity of the default-mode (and other resting state networks) occurs primarily in the low frequencies (<0.1Hz; (Cordes et al., 2001)), so the spatial extent of the network may be insensitive to small latency differences. Indeed, the autocorrelation of the PCC ROI timeseries has a width of several seconds (Figure 12). We expect that when examining correlations with a reference timeseries having substantial high-frequency components, latency correction may have a higher impact on the resulting networks. Indeed, the impact of latency correction on activation maps in our event-related WM task was greater. Post-correction Granger causality analysis in one subject also revealed significant causal connections across the working memory network that were absent prior to the correction. While anecdotal, this result illustrates the potential impact of vascular latency differences on quantifying causal relationships between regions.
Quantifying relative hemodynamic latency using BH task may contribute to a wide range of fMRI modeling studies. For example, in analyses of causality where it may be difficult or impossible to use load-dependent modulation (e.g. resting state), a causality mapping applied to time series in the BH task may serve as a “control” analysis, indicating whether an observed causal relationship in the task of interest could in fact be due to vascular delay (e.g. Figure 10b). This is particularly critical for connectivity methods that do not include an explicit forward model, such as Granger causality and structural equation modeling. Latencies measured with a BH task might also be entered as parameters in deconvolution methods and biophysical models (e.g. dynamic causal modeling (Friston et al., 2003)). DCM addresses the hemodynamic confound in neural connectivity modeling by allowing specification of a hemodynamic forward model, and simultaneously estimating both connectivity and region-specific HRF paramters (including latency) (Friston et al., 2003; Stephan et al., 2007). Our approach may thus be viewed as complementary, or confirmatory, to such methods. Importantly, our approach provides a means of gauging relative hemodynamic delays across regions without relying on a detailed (and suitable) specification of model parameters; indeed, unless one is certain of the model, data-driven approaches can yield misleading results, and independent measurements of latency may prove useful before DCM, ICA, and other non-hypothesis driven methods are employed.
Finally, we note that with the exception of cerebellar areas and major vessels, relative latency differences across the brain tended to be small, and the correction had only a minor impact on the resting state networks and event-related tasks examined in the current study. Therefore, although vascular latency correction is critical for studies in which the precise timing of regional responses is of interest (e.g. causality analysis), for many applications it may suffice to employ a sensory task for inferring hemodynamic latency in most other non-sensory regions.
This research was supported by NIH grants P41-RR09784 to GG and T32-GM063495 to CC.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.