Simultaneous EEG-fMRI Reveals Temporal Evolution of Coupling between Supramodal Cortical Attention Networks and the Brainstem

Behavioral paradigm.

Seventeen subjects (six females; mean of 27.7 years; range, 20–40 years) participated in three runs each of analogous visual and auditory oddball paradigms. The 375 (125 per run) total stimuli per task were presented for 200 ms each with a 2–3 s uniformly distributed variable intertrial interval and target probability 0.2. The first two stimuli of each run were constrained to be standards. For the visual task, the target and standard stimuli were, respectively, a large red circle and a small green circle on isoluminant gray backgrounds (3.45° and 1.15° visual angles). For the auditory task, the standard stimulus was a 390 Hz pure tone, which was selected to lie within a trough of the scanner sound frequency spectrum, and the target sound was broadband. This broadband “laser gun” sound was chosen such that EEG discriminator performance matched that of the visual task. Because our study focused on task-related attentional states, subjects were asked to respond to target stimuli, using a button press with the right index finger on an MR-compatible button response pad. Stimuli were presented to subjects using E-Prime software (Psychology Software Tools) and a VisuaStim Digital System (Resonance Technology), comprising headphones and 600 × 800 goggle display. All subjects gave informed consent following the protocol of the Columbia University Institutional Review Board.

Simultaneous EEG and fMRI data acquisition.

A 3 T Philips Achieva MRI scanner (Philips Medical Systems) was used to collect functional echo-planar image (EPI) data with 3 mm in-plane resolution and 4 mm slice thickness. We covered the entire brain by obtaining 32 slices of 64 × 64 voxels using a 2000 ms repetition time and 25 ms echo time. We also acquired a single-volume high-resolution (2 × 2 × 2 mm) EPI image and a 1 × 1 × 1 mm spoiled gradient recalled image for each subject for purposes of registration.

We simultaneously and continuously recorded EEG using a custom-built MR-compatible EEG system (Goldman et al., 2009; Sajda et al., 2010; Walz et al., 2013), with differential amplifier and bipolar EEG cap. The caps were configured with 36 Ag/AgCl electrodes including left and right mastoids, arranged as 43 bipolar pairs. Bipolar pair leads were twisted to minimize inductive pickup from the magnetic gradient pulses and subject head motion in the main magnetic field. This oversampling of electrodes ensured data from a complete set of electrodes even in instances when discarding noisy channels was necessary. To enable removal of gradient artifacts in our offline preprocessing, we synchronized the 1-kHz-sampled EEG with the scanner clock by sending a transistor–transistor logic pulse to a field-programmable gate array card (National Instruments) at the start of each of 170 functional image acquisitions. All electrode impedances were kept below 20 kΩ, including 10 kΩ resistors built into each electrode for subject safety. A comprehensive description of the hardware, along with the preprocessing and analysis methods described throughout the remainder of Materials and Methods, can be found in the book chapter by Sajda et al. (2010).

EEG data preprocessing.

We performed standard EEG preprocessing offline using MATLAB (MathWorks) with the following digital Butterworth filters: 1 Hz high pass to remove direct current drift, 60 and 120 Hz notches to remove electrical line noise and its first harmonic, and 100 Hz low pass to remove high-frequency artifacts not associated with neurophysiological processes. These filters were applied together in the form of a linear-phase finite impulse response filter to avoid distortions caused by phase delays.

Ballistocardiogram (BCG) artifacts are more challenging to remove, because they share frequency content with EEG activity. Existing BCG removal algorithms cause loss of signal power in the underlying EEG, so we performed single-trial classification (described in Materials and Methods, Single-trial analysis of EEG) on the data before BCG artifact removal. This is a justified practice because our classifier identifies discriminating components that are likely to be orthogonal to BCG. To compute scalp topographies of these discriminating components, BCG artifacts were removed from the continuous gradient-free data using a principal components analysis method (Goldman et al., 2009; Sajda et al., 2010). First, the data were low passed at 4 Hz to extract the signal within the frequency range in which BCG artifacts are observed and then the first two principal components were determined. The channel weightings corresponding to those components were projected onto the broadband data and subtracted out. These BCG-free data were then re-referenced from the 43 bipolar channels to the 34-electrode space to calculate scalp topographies of EEG discriminating components (described in Materials and Methods, Single-trial analysis of EEG).

We extracted 1000 ms stimulus-locked epochs, with baseline removal on the 200 ms before stimulus, from both the 43-channel datasets and the 34-electrode re-referenced dataset. By visual inspection, we discarded trials containing motion or blink artifacts, evidenced by sudden high-amplitude deflections, and also those with incorrect responses, identically for both datasets. This left ∼95% of the target trials remaining.

Single-trial analysis of EEG.

Independently for the auditory and visual data, we applied a linear classifier to the 43-channel EEG signal amplitude using the sliding window method of Parra et al. (2005). Specifically, we found a projection of the multidimensional EEG signal, x_i(t), where i = {1… T} and T is the total number of trials, within a short time window that achieves maximal discrimination between target and standard trials. All time windows had a width of N = 50 ms, and the window center, τ, was shifted from 0 to 1000 ms relative to stimulus onset, in 25 ms increments. We used logistic regression to learn the 43-channel spatial weighting, w(τ), that maximally discriminated conditions, arriving at the projection, y_i(τ), for each trial i and a given window τ. Note that we use the average projection in each temporal window, τ, which tends to be less sensitive to noise (Parra et al., 2005). The result is a hyperplane, for each τ, that maximally discriminates the target versus standard trials given the EEG amplitudes for each trial (Fig. 1A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Model design. The entire method was applied independently for multiple temporal windows, τ, of EEG data spanning the trial. Shown is an example for a given τ. A, We first estimate w, which is a linear weighting on the EEG sensors that maximally discriminates stimulus conditions: targets (red) vs standards (green). This determines a task-relevant projection of the data, in which the distance to the decision boundary reflects the decision certainty of the classifier. B, From w, we compute ȳ_i, which is the demeaned classifier output, ȳ_i = w^T x_i − y_mean for each target trial, i. We treat ȳ_i as an index of attentional state on each target trial. Note that ȳ_i contains no discriminating information about target versus standard trials and instead represents the STV of the task-relevant EEG projection around the mean of all targets. C, Given the ȳ_i values and their corresponding stimulus-onset time points, we build regressors for an fMRI GLM, including event-related average BOLD response to targets (D), RT variability (E), and residual EEG variability (F). Similar regressors were included for the standard stimuli, along with motion parameters and temporal derivatives. All were convolved with the canonical HRF.

The classifier output, y_i(τ), represents the confidence of the classifier in its prediction based on the training data and the model. We assessed classifier performance with the area under the receiver operating characteristic curve (Green and Swets, 1966), denoted AUC, using leave-one-out (LOO) cross-validation (Duda et al., 2000). AUC was calculated for multiple temporal windows, enabling observation of the temporal progression of task-relevant components and localization of the event-locked time with maximal discrimination between conditions. To obtain a significance threshold for the AUC values, we used a permutation test in which we randomly permuted the trial labels and ran the classifier using the LOO approach. We repeated this procedure 1000 times for each subject to generate a distribution of AUC values from which we computed the null hypothesis distribution for AUC values and the corresponding AUC threshold for p = 0.01. For each window τ, we also generated the forward model a(τ), where y_i(τ) is now organized as a vector y(τ), where each row is from trial i, and x_i(t) is organized as a matrix, X(τ), where rows are channels and columns are trials, all for time window τ. These forward models can be viewed as scalp plots and interpreted as the coupling between the discriminating components and the observed EEG (Parra et al., 2005; Goldman et al., 2009; Sajda et al., 2010; Walz et al., 2013). Because we estimate a linear weighting for temporal windows spanning the trial, we are able to track the progression of the spatial components across time.

fMRI data preprocessing.

Using FSL (FMRIB Software Library; Smith et al., 2004), we performed bias-field correction on all images to adjust for field distortion artifacts caused by the EEG wires. We then performed slice-timing correction, motion correction, 0.01 Hz high-pass filtering, and 5 mm full-width half-maximum spatial smoothing on the functional data. Motion correction provided motion parameters, which were later included as confounds in the GLM. Functional and structural images were registered to a standard Montreal Neurological Institute (MNI) brain template after brain extraction, and each subject's image registration was checked manually to ensure proper alignment.

Traditional fMRI analysis.

We first ran a traditional fMRI analysis, using event-related and RT variability regressors in our GLM. The event-related regressors comprised boxcar functions with unit amplitude and onset and offset matching that of the stimuli. RT variability was modeled using unit amplitude boxcars with onset at stimulus time and offset at RT, and these were orthogonalized to the event-related regressors. Orthogonalization was implemented using FSL, which uses the Gram–Schmidt procedure (Strang, 2003) to decorrelate the RT regressor from all other event-related regressors. All regressors were convolved with the canonical hemodynamic response function (HRF), and temporal derivatives were included as confounds of no interest. An event-related targets versus standards contrast was also constructed. A fixed-effects model was used to model activations across runs, and a mixed-effects approach was used to compute the contrasts across subjects.

We combined individual runs of the experiment at the subject level in three ways: (1) including only auditory runs; (2) including only visual runs; and (3) including all six auditory and visual runs together. Each of these three analyses was carried through to the group level. Statistical image results for these traditional analyses were thresholded at z > 2.3, and clusters were multiple-comparison corrected at p < 0.05 (Worsley, 2002). Because it is possible for significant correlations to be driven much more strongly by a subset of the data (in this case by one task), we used the subject-level auditory-only and visual-only results to compute paired t tests and thus identify any brain regions with BOLD correlates that differed significantly across auditory and visual runs. We discarded any voxels that achieved an uncorrected voxelwise p < 0.01 from the statistical maps of the combined analysis to ensure that our reported findings are common to both auditory and visual tasks.

Because pulsation artifacts in the ventricles have been attributed to false activations near the BS (Astafiev et al., 2010), we similarly excluded any voxels that were correlated with the BCG pulsation artifact present in the EEG. These voxels were determined with a p < 0.01 uncorrected threshold, after incorporating the pulse timing into a GLM and carrying through to group level. The pulse timing was obtained using a peak detection algorithm on the EEG signals at temporal electrodes (T7, T8, CP5, CP6), which are the sites most strongly affected by such artifacts.

EEG-based fMRI analysis.

For the single-trial variability (STV) fMRI analysis, we modeled the variability of the neural response using an additional two regressors: one each for targets and standards. These EEG-based regressors were designed with 100 ms duration, centered on the classifier training window. The STV regressor height was modulated using the de-meaned output, ȳ_i = y_i − y_mean, of the EEG discriminator for each trial, i (Fig. 1). These regressors were convolved with the HRF and orthogonalized with respect to all traditional regressors, with temporal derivatives included as confounds. It was especially important to regress out the RT variability, because RT is known to be negatively correlated with attention and our aim was to study variability in task-related attention that cannot be detected using an external measure. This entire analysis was run independently for all EEG training windows exceeding a mean AUC value of 0.75 for both tasks, which is a common psychophysical threshold used in signal detection theory and highly conservative with regard to EEG discrimination. We focused on within-class variability, using only the task-relevant (target) stimuli STV statistical maps in our results interpretation.

We used the following randomization method to generate a null distribution and correct for multiple comparisons. We ran the full STV fMRI analysis described above after permuting the y_τ = 450 ms values randomly within each class. One hundred such permutations were done for each subject, from which we generated 10,000 group-level statistical maps via random samplings from the subject-level results. Stelzer et al. (2013) demonstrated the validity of this oversampling method. This null distribution was used to generate familywise error (FWE)-corrected positive and negative activation p value maps within the BS (one map for each of the EEG windows).

Cluster-size thresholding is commonly used to correct for multiple comparisons in cortical regions, because true activations are assumed to be widespread and FWE correction is often too strict for such a large number of voxels. To avoid an arbitrarily selected cluster-size threshold in our cortical analysis, we applied threshold-free cluster enhancement (TFCE) to account for variable z-score magnitude within each cluster (Smith and Nichols, 2009). After also applying TFCE to all maps comprising the null distribution, we used FWE correction within a gray matter mask. All corrected maps were thresholded at p = 0.05.

As in the traditional analysis, we ran the EEG STV fMRI analysis separately three ways: (1) only auditory runs; (2) only visual runs; and (3) with auditory and visual combined at the subject level. The latter method was of primary interest, because it was used to investigate BOLD correlates that generalize across both sensory modalities. As described for the traditional fMRI analysis (see below, Traditional fMRI analysis), we used the auditory-only and visual-only results to determine voxels with activation that significantly differed between the two tasks. These were then excluded from the statistical maps of the combined analysis, as were any voxels that correlated with BCG pulse timing. No more than 5% of the initial voxels were eliminated. This ensured that our reported findings generalize across the auditory and visual domains.

Dynamic causal modeling of fMRI data.

The EEG-based fMRI analysis described above provided timing information about the BOLD correlates of endogenous attentional modulations, but it did not allow inference of directed relationships between these regional activations. We investigated interactions between our temporally specific EEG variability correlates by performing an effective connectivity fMRI analysis. In particular, we implemented single-state linear dynamic causal modeling (DCM; Friston et al., 2003; Stephan et al., 2010) using DCM10 in SPM8.

We used the results of our EEG-based fMRI analysis (see Fig. 6) to select five regions of interest (ROIs); these were chosen based on suprathreshold significance and particular interest in relation to previous studies of attentional networks: (1) BS; (2) ACC; (3) right ACC (rACC); (4) right middle frontal gyrus (rMFG); and (5) right orbitofrontal cortex (rOFC).

The activations were spread across four different EEG time windows (i.e., correlated with variability in temporally specific EEG discriminating components), and some differed in sign of correlation so their time series were not necessarily intrinsically correlated. Each ROI was a sphere with 4 mm radius, centered on the peak voxel of the group-level EEG-based GLM results. Time series were extracted from individual subjects' preprocessed functional data in MNI space (downsampled to 3 × 3 × 3 mm for computational efficiency) by estimation of the first principal component within each ROI. Auditory and visual runs of the experiment were combined because the focus of our study was task-related attentional modulations common across sensory input modalities.

We designed 12 related models (see Fig. 8) to investigate intrinsic directed connectivity among the five ROIs, focusing on cortical projections to and from the BS. In Models 1–9, we varied the direction of the ACC–BS connection (including one-way and bidirectional connections) and the stimulus input region (BS, ACC, or both). Directionality of the other connections was based on known structural projections in primates (Aston-Jones and Cohen, 2005), including a directed connection from the BS to rMFG and bidirectional connections between BS and rACC, as well as BS and rOFC. Models 10–12 used timing information from the EEG-based analysis in specification of all the connections, only allowing connections to be directed forward in time. Input region was varied similarly to the other models.

We used fixed-effects Bayesian model selection (BMS) to compare these 12 models both on a single-subject level and at the group level. BMS balances model fit and complexity, thereby selecting the most generalizable model. It estimates the relative model evidence and provides a distribution of posterior probabilities for all of the models considered. We used fixed effects because focused attention to a simple task involves basic physiological processes that are expected to be common to a healthy subject population (Stephan et al., 2009, 2010). We also compared families of similar models (Penny et al., 2010); the model space was divided into three families based on the direction of connectivity between the BS and ACC. Finally, we used Bayesian parameter averaging (BPA) to provide group-mean intrinsic connection strengths and their probabilities for the winning model.