Functional Source Separation for EEG-fMRI Fusion: Application to Steady-State Visual Evoked Potentials

1. Introduction

An organism's survival depends on its ability to quickly and adaptively respond to environmental challenges and opportunities. To accomplish this task, the mammalian brain detects and stores the predictive value of recurring environmental signals with respect to dangerous or rewarding outcomes (Sokolov, 1963). Established associative networks that link specific stimuli to representations of biological significance and motor action are the mutual interest to Neuroscientists and Neurorobotists (Falotico et al., 2017; Oess et al., 2017). An extensive literature has demonstrated that the neural and hemodynamic amplification of threat, relative to neutral, cues is paralleled by a host of behavioral effects such as facilitated detection (Öhman and Soares, 1998), identification (Anderson, 2005), and greater perceptual vividness (Todd and Thompson, 2015).

Hemodynamic and electrophysiological time series represent the two major measurement modalities used in cognitive neuroscience research today. They convey complementary, unique information, with each modality possessing sensitivity to different facets of brain physiology. Above and beyond their well-established differences in temporal and spatial resolution, fMRI-BOLD signals, and electrical field potentials often display divergent responses to experimental manipulations, even in situations where they are considered indices of the same large-scale neural process. Thus, characterizing the relationship between hemodynamic and electrophysiological time series has important implications for neurophysiological studies of the human brain. Consistent with recent empirical studies into the nature of the BOLD signal (Logothetis et al., 2001; Logothetis, 2015), and with the biophysics of EEG (Nunez and Srinivasan, 2006), the EEG-fMRI fusion approach examined in the present study assumes that co-variation between EEG and fMRI signals accounts only for fraction of the variance of each measure (Herrmann and Debener, 2008). To facilitate cross-validation and precise prediction of expected patterns of co-variation, we use non-invasive simultaneous EEG-fMRI recordings during a fear conditioning task in which strong and sustained visuocortical responses are elicited. Specifically, we combine fMRI with steady-state visual evoked potentials (ssVEPs), elicited by stimuli that are rapidly and periodically modulated in luminance or contrast. Using standard psychophysiological techniques such as EEG, these responsing signals can be recorded at the surface of the scalp as an oscillatory waveform that has the same fundamental frequency as the driving stimulus (Regan, 1989). Readily quantified at the level of single trials and of known origin in extended visual cortical areas, this approach allows us to capture different aspects of co-variation in simultaneously EEG and fMRI (Liu et al., 2012): the correlation may reflect responses to experimental events seen in both modalities, mixing trivial effects of joint reactivity of EEG and fMRI to salient external events with correlations of interest such as condition-specific co-variation, or co-variation with various cognitive states. In addition, and of particular relevance in the absence of experimental stimulation, specific portions of the co-variance between modalities may reflect neural-hemodynamic dependencies in spontaneous, ongoing brain activity. In previous work, we have described a method for quantifying the spatial and temporal correlations between the time series of fMRI BOLD and suitable indices derived from the EEG signal (Ji et al., 2018). The present work examines in more detail the benefits of ssVEP preprocessing using Functional Source Separation (FSS) as an optimization step.

The ssVEP possesses a high signal-to-noise ratio in addition to other desirable qualities for EEG-fMRI fusion: The frequency of the ssVEP is precisely known and can therefore be reliably separated from noise and quantified in the frequency domain (Spekreijse et al., 1985; Regan, 1989; Wang et al., 2007). The high signal to noise ratio of ssVEPs allows the robust quantification of visual neuronal population responses at the level of individual trials, of particular interest in the context of EEG-fMRI studies. In terms of neurocognitive function, fluctuations in ssVEP amplitude are sensitive to stimulus content (e.g., Keil et al., 2003) as well as task instructions (e.g., Müller et al., 1998; Attar et al., 2010). Generators of the ssVEP are in the extended visual cortex (including occipital, temporal, and parietal cortices; Müller et al., 1997), with most pronounced contributions from retinotopic areas (Wieser and Keil, 2011), but also from cortices higher in the visual hierarchy, along the ventral and dorsal streams (Di Russo et al., 2007). The ssVEP technique thus represents a robust and reliable method for non-invasively isolating population-level neuronal responses at low levels of the traditional visual hierarchy (Spekreijse et al., 1985; Regan, 1989), very well-suited for joint validation of EEG-fMRI analyses, in which it has been previously used (Sammer et al., 2005). Existing pipeline for ssVEP-fMRI analyses have shown high variability across participants and surprisingly small co-variation of ssVEP amplitude and BOLD magnitude in visual areas. Given that flickering or contrast-reversing stimuli evoke very strong BOLD and electrophysiological responses, it is conceivable that previous findings reflect noise in the single-trial estimation of ssVEPs rather than a lack of co-variation. The present study illustrates the application of a new preprocessing pipeline for ssVEP-fMRI fusion and presents its outcomes regarding group-level findings, including the correspondence of ssVEP-fMRI covariation maps among individual participants (see Figure 1).

FIGURE 1

Figure 1. Preprocessing pipeline of ssVEP-fMRI fusion using Functional Source Separation (FSS) as an optimization step, timing information which enhances the identification and separation of relevant signals in electrophysiological time series while maximizing negentropy constraint in the meantime, the summarized ssVEP time series then go for the spatial-temporal correlation analysis with a proper temporal delay.

EEG quality is a crucial issue when conducting joint EEG-fMRI analyses. Blind source separation (BSS) represents a group of techniques for extracting signals in presence of noise, by exploiting the statistical properties of sources. For example, independent component analysis (ICA) recovers independent sources that contribute to form the measured signal, under the assumptions that they are linearly mixed, and measured with an array of spatially distributed sensors (Comon, 1994; Bell and Sejnowski, 1995; James and Hesse, 2004). ICA is widely used in neurophysiological time series analysis (Makeig et al., 2002; Delorme and Makeig, 2004) and is increasingly the standard method when performing single-trial based EEG-fMRI fusion (Debener et al., 2005; Eichele et al., 2005; Bénar et al., 2007; Mulert et al., 2008) following removal of gross MRI artifacts. It has also been successfully used for removal of eye blinks, eye movements, and electrode artifacts (Debener et al., 2006; Mantini et al., 2007), denoising muscle artifacts (Chen et al., 2017, 2018) and MR noises (Acharjee et al., 2015). Several novel techniques have been proposed for multiset and multimodal fusion analysis (Chen et al., 2013, 2016). Instead of an entirely blind searching approach that relies purely on statistical properties of the signals, the current study used Functional Source Separation (FSS). This method uses timing information to constrain the solutions, that satisfy physiological assumptions. FSS thus represents an optimization step in EEG-fMRI fusion, with the goal of improving the signal quality of EEG data recorded in the MRI environment. FSS is a semi-blind extension of ICA that incorporates prior knowledge about the responses of interest (Barbati et al., 2006). Thus, FSS enhances the identification and separation of relevant signals in electrophysiological time series while aiming to maximize the neurophysiological validity that is typically associated with using the averaged signal (Ostwald et al., 2010, 2011; Porcaro et al., 2010, 2011; Naik and Wang, 2014). In the present pipeline, following removal of gradient and cardioballsitic artifacts with standard techniques, data were further processed using FSS.

The present paper illustrates the usage of this FSS based approach for quantifying ssVEP-fMRI co-variation. Using a fusion algorithm that we have previously described (Ji et al., 2018), we applied the FSS preprocessing step to data from a 3-phase aversive conditioning paradigm that included a habituation, acquisition, and extinction phase. EEG and fMRI were recorded simultaneously while participants viewed periodically (10Hz) and rapidly phase-reversing Gabor patches (sine-wave gratings), evoking steady-state visual potentials (ssVEPs).

2. Materials and Methods

2.1. Participants

Participants consisted of 10 undergraduate and graduate students (4 male, average age 20.8 years) who, after giving written informed consent, participated on a volunteer basis or received course credit in the General Psychology course taught at the University of Florida. All participants were screened for metallic implants, claustrophobia, and history of seizure episodes. Female participants self-administered a pregnancy test prior to participation. Data of one additional subject was excluded from the analysis for poor quality of EEG data, defined as retaining <50% of the trials after artifact screening. All procedures were approved by the institutional review board of the University of Florida and were consistent with the Declaration of Helsinki on studies with human participants.

2.2. Stimuli and Experiment Setup

Stimuli were Gabor patches, i.e., sinusoidal gratings multiplied with a Gaussian envelope. They were oriented at either 15 or 345° relative to the vertical meridian. Gabor patches were shown at a maximum Michelson contrast of 95% (maximum = 110cd/m²; minimum = 2.1cd/m²) on a dark gray background, at a spatial frequency of 0.45 cycles per degree, and a horizontal visual angle of 15.5°. Participants viewed all stimuli on a mirror placed on the MR head-coil positioned 8.5 cm from eyes, allowing them to see an MR-compatible monitor that was placed outside the scanner bore. The stimuli reversed their phase every 100 ms to evoke a ssVEP having its second harmonic at 10 Hz. Consistent with previous studies of phase-reversal ssVEP (see Keil et al., 2013; Norcia et al., 2015 for a review), we analyzed the second harmonic response, i.e., 10 Hz.

The data presented in this report were recorded from a differential aversive conditioning study in which Gabors of one orientation (but not the Gabor patches with the other orientation) were occasionally paired with an electric shock (see Petro et al., 2017, for details). For the habituation block, participants were instructed that they would not feel any shock but to fixate on the patterns. During the acquisition block, participants were informed that they would intermittently feel a cutaneous electric shock during the experiment but were not instructed as to the contingencies of the shock administration. The extinction phase was also uninstructed, such that participants were not told that no more shocks were to be given. Each participant was instructed to remain still while in the scanner and to maintain fixation on the center of the screen. All stimulus presentation was controlled using E-Prime software (Psychology Software Tools, Pittsburgh, PA). Cutaneous shocks were administered using an STMISO Stimulation Isolation Adapter (BIOPAC Systems, Inc., Goleta, CA) with MRI-compatible skin electrodes.

The data reported here include 40 total trials per phase per participant. Each trial consisted of one of the two gratings being presented for 5,100 ms, during which its phase was alternated every 100 ms. An inter-trial interval (ITI) consisted of an initial gray cross (37.5cd/m²; 1° of visual angle) presented in the middle of the screen for a random duration between 0 and 8 s followed by a white cross (149.0 cd/m²) for a duration of 3 s, immediately preceding trial onset with Gabor patch presentation.

EEG data were recorded on a 32-channel MR compatible system (BrainAmp MR plus amplifier by BrainProducts). This system consisted of thirty-one Ag/AgCl electrodes placed on the head according to the 10–20 system and one cardioballistic (CB) electrode placed on the upper back to record heart-beats. The reference was positioned at FCz, the ground electrode was placed 1 cm anterior to Oz. Impedances were reduced to below 10 kΩ for all scalp electrodes and below 50 kΩ for the cardio electrode, as suggested by the Brain Products manual. According to the Brain Vision user manual, a sampling rate of 5 kHz EEG permits the TDD (Template Drift Detection) algorithm the highly accurate detection of drifts and shifts in typical MR data set, therefore EEG data were recorded online at 5 kHz and digitized to 16 bit, while the digitized data is transferred via a fiber-optic cable to the recording computer. The system was synchronized to the internal clock of the scanner, and synchronized event marker has been added for further signal processing.

MRI data were collected with a 3T Philips Achieva scanner. An Avotec Silent Scan headphone system was used to diminish gradient noise. Data were acquired during gradient-echo echo-planar imaging sequence [echo time (TE), 30 ms; repetition Time (TR), 1.98 s; flip angle, 80°; slice number, 36; field of view, 224 mm; voxel size, 3.5 × 3.5 × 3.5mm³; matrix size 64 × 64]. The first four functional scans were discarded to allow for scanner stabilization. Slices were acquired in ascending order oriented parallel to the plane connecting the anterior and posterior commissure during an 1,850 ms interval with 130 ms between each TR, during which no images were collected and allowed visual inspection of the EEG data during recording when the MR gradient artifact was absent. A T1-weighted high-resolution structural image was obtained after completion of all functional scans.

2.3. Artifact Handling and Preprocessing

Brain Vision Analyzer 2.0 (Brain Products) Software was used to remove artifacts caused by the magnetic gradients and heart-beats. Gradient artifacts are due to the switching of magnetic resonance (MR) gradients that is necessary to collect MRI data. The ballistocardiogram (BCG) artifact is related to the pulsatile movement of blood and the pulsatile movement of electrodes adjacent to large blood vessels. These two sources of noise are relatively independent. The removal of magnetic gradient artifacts followed a modified version of the algorithm developed in Allen et al. (2000), in which an artifact template constructed from a sliding window of 41 consecutive volumes is subtracted from the EEG data according to the event markers sent from the scanner's internal clock. Cardioballistic artifacts were removed with a standard algorithm as well (Allen et al., 1998). This method detected R peaks in the CB electrode and created a template from 21 consecutive heartbeats, which were subtracted from the continuous EEG data. Particular attention was given to reliable ECG peak detection: Misaligned and falsely detected R peak markers were corrected manually according to the peak of the correlation curve. Data was then down-sampled to 250 Hz. To isolate non-cerebral signals such as eye blinks, ICA was applied, and components with distinct peaks at sensors FP1 and FP2 removed (Delorme and Makeig, 2004). Then, the data were referenced to a common average, removing any global, non-cerebral signals that were common to all sensors.

Preprocessing of BOLD fMRI data was completed using SPM12 (Eickhoff et al., 2005, freely available at http://www.fil.ion.ucl.ac.uk/spm/). We followed the standard preprocessing routines suggested by SPM: Time differences were compensated by slice timing correction. Head movements were estimated by realigning each scan to match one representative scan with rigid transformation. Images were normalized and registered to the Montreal Neurological Institute (MNI) space. Functional volume images were re-sampled to a spatial resolution of 3 × 3 × 3 mm³. Images were smoothed using a Gaussian kernel with a full-width at half-maximum of 6 mm. Low-frequency temporal drifts were removed from the BOLD data using a 1/128 Hz high-pass filter.

2.4. FSS With Prior Information

Electroencephalography (EEG) reflects postsynaptic activity in large populations of cortical neurons that have similar spatial orientation (Nunez and Srinivasan, 2006). EEG signals are often assumed to be generated by a linear mixing process of many sources of activity (i.e., many such populations), which may overlap spatially and temporally and are distorted due to the effect of volume conduction. For the ssVEP, the known origin of the signal is located in lower-tier (calcarine fissure) and extended visual cortical areas. Because of the desirable features of the ssVEP (known frequency, high signal-to-noise ratio) robust ssVEP signals can be obtained in the fMRI scanner environment. As illustrated by Figure 2, the EEG epochs were extracted in the range of 1 s pre-stimulus and 6 s post-stimulus onset. The resulting data segments of 7 s were transformed into the frequency domain (top left). In the resulting spectrum, the maximum absolute value of the Fourier coefficients was observed at the driving frequency of 10 Hz, most pronounced at posterior electrodes (topographic map, top right). Furthermore, the Grand Mean time domain average across all trials and subjects shows enhanced evoked oscillatory activity (bottom) during the period when the phase-reversing stimulus was on, compared to baseline. FSS aims to maximize the neurophysiological validity that is associated with this latter property of the averaged signal. Therefore, we use timing information to constrain the solutions, while also considering the independence source constraint, which enhances the identification and separation of relevant signals in electrophysiological time series. The modified contrast function is:

$\begin{array}{ll}F(\text{w})=J(\text{w})+\lambda H(\text{w})& (1)\end{array}$

The FSS optimization tends to find a projection vector w that maximizes the above contrast function. J can be any function for ICA. Here, we adopted the criteria of fastICA (Hyvarinen, 1999), aiming to minimize the Gaussianity of the results (equivalently, to find the separation vector in which the negentropy is maximized), since the observed mixed signals will tend to have more Gaussian amplitude distributions. The function H accounts for prior information we have on sources, which give preference to solutions that tend to satisfy physiological assumptions. Parameter λ is used to weigh the two parts of the contrast function, chosen to be 10 for ssVEP recovery.

FIGURE 2

Figure 2. Grand Mean time series of filtered (0.5–30 Hz) EEG time-series (Oz) time locked to trial onset and averaged across all trials and subjects (bottom). Oscillatory activity driven by stimuli picture evoked flicker stream (10 Hz) are the dominant signals as illustrated by the frequency domain spectrum (top left), most strongly represented by posterior electrodes in the topographic map (top right).

The function J serves as a generalization of the cumulant-based approximation of negentropy

$\begin{array}{ll}J(\text{w})={\left[\text{E}\left\{\text{G}({\text{w}}^{\text{T}}\text{x})\right\}-\text{E}\left\{\text{G}(\text{v})\right\}\right]}^2& (2)\end{array}$

where x is the pre-processed EEG and v is a standardized Gaussian variable. The variance of w^Tx must here be constrained to unity. The following choices of G have proven useful according to Hyvarinen (1999), Hyvärinen and Oja (2000)

$\begin{array}{ll}G(u)=\frac{1}{a_1}log(cosh(a_{1}u))& (3)\end{array}$

where 1 ≤ a₁ ≤ 2 is some suitable constant, we here chose a₁ = 1.

To facilitate the FSS optimization, a zero-phase narrow bandpass filter centered at the target frequency 10 Hz with plus and minus 0.5 Hz was applied. The stimulus response is defined as

$\begin{array}{l}{R}_{stim}=\frac{1}{t_e-t_s}{\displaystyle {\int }_{t_s}^{t_e}}abs(Hilbert(EA(t)))dt\\ -\frac{1}{0-t_b}{\displaystyle {\int }_{t_b}^0}abs(Hilbert(EA(t)))dt& (4)\end{array}$

with the evoked activity, EA, computed by averaging signal epochs of the estimated source w^Tx, which lasted for 5 s for the post-stimulus reactivity and 1 s for the pre-stimulus baseline. Accordingly, the upper and lower limits of the integral are set as t_s = 0, t_e = 5 and t_b = −1. The operation abs(Hilbert(·)) extracts the envelope of the evoked activity (Bracewell and Bracewell, 1986; Huang, 2014). Following the former work (Barbati et al., 2006) which apply FSS on MEG data, a saturation function is adopted in function H and tuned to suit the ssVEP study:

$\begin{array}{ll}H=\psi (R_{stim},k)& (5)\end{array}$

where

$\psi (R_{stim},k)=\{\begin{array}{ll}{R}_{stim}/k\hfill & \text{when\hspace{0.17em}}{R}_{stim}\le \text{k}\hfill \\ 1\hfill & else\hfill \end{array}$

where k is the parameter measuring the required minimum response, so as to define an admissible region where the optimization is only driven by J when response is greater than k. For the present application, the maximum R_stim for each phase is computed separately with fastICA, and assigned to k. Since prior information about the sources can only be described by a non-differentiable function, simulated annealing is employed for global optimization (Kirkpatrick et al., 1983; Ingber, 2000). To ensure the stability and reliability of the optimization procedure, the contrast function for all three phases were added up and jointly optimized.

The approach described above gives us the possibility to extract only one component that maximizes the functional behavior in agreement with the functional constraint.

2.5. Quantifying the SNR of Extracted ssVEP

In order to evaluate the performance of FSS optimization, the signal-to-noise ratio is determined by comparing the signal level during the stimuli on period with the signal level without stimulus. The pre-stimulus baseline from −1 to 0 s was set as the time interval of background/baseline noise, during which the participant was instructed to maintain fixation on a white cross, preceding the onset of the Gabor patch. In addition, a period from the inter-trial baseline (5.6–7 s after stimulus onset) was also included as a manipulation check, quantifying the extent to which the extracted ssVEP oscillation attenuates after stimulus offset. Since this part of background information is blind to FSS optimization, it represents an intuitive test for neurophysiological validity.

2.6. ssVEP-BOLD Fusion

To assess the spatial-temporal correlation between EEG and fMRI responses recorded simultaneously, the recovered ssVEP source was averaged across time segments of each fMRI scan, representing the intensity of neuronal activities, resulting in a time series with the same temporal resolution of the BOLD. As in the procedure of FSS optimization described in methods, the proper scaling (including sign) of the source signals is non-recoverable, we are mostly interested in the predictability reflects the proportion of the variance in the BOLD signal that is linearly predicted from ssVEP time series. Following the previous work (Wang and Zheng, 2014; Ji et al., 2018), we use the bivariate case of cross multivariate correlation coefficient (xMCC) for correlation analysis between EEG and fMRI, which is equal to the absolute value of the cross correlation coefficient. The metabolic (BOLD) events are slower and temporally lagged to the neuro-electric activity (EEG). Based on the large body of research on stimulus-BOLD latency performed in the visual cortex (Buxton et al., 1998; Serences, 2004; Penny et al., 2011), we adopt the consensus setting of 4 s as the latency of BOLD relative to neural events.

2.7. Statistical Inference: Permutation Testing

To assess the statistical significance of the correlation results, we applied permutation tests to determining the threshold for rejecting the null hypothesis of no linear relation based on shuffled data, and assumed statistical significance (permutation controlled in brain regions) if p-value is less than 0.05 for single participant. Given the small size of the present sample, the robustness of a given ssVEP-fMRI correlation across participants was quantified by measuring inter-subject overlap of voxels showing significant correlations after thresholding by a permutation test, using a binomial test. To this end, we calculated the probability of overlap at a given voxel under the null hypothesis that significant correlations occurred in different voxels and different participants. The resulting p-values were then corrected by the number of voxels considered overall, resulting in a threshold of 0.5 at a given voxel (i.e., five out of 10 participants showing permutation controlled, significant correlations at that voxel). These results were compared to the average correlation map, to establish the validity of this conservative statistical approach.

3. Results

3.1. EEG-ssVEP Source Extraction

The SNR index was evaluated for each participant and experimental phase separately, and compare with ICA results optimized for each phase separately. As a referencing index, in Table 1 we reported the spectrum signal-to-noise ratio (SNR_RAW), which divides the power of ssVEP at 10 Hz by the mean power of neighboring bins in the spectral range between 0.5 and 30 Hz, together with ssVEP amplitude (ssVEP_AMP) (Keil et al., 2008), by averaging the EEG traces in the moving windows in the time domain, which was then submitted to a Discrete Fourier Transform to get the 10 Hz ssVEP amplitude.

TABLE 1

Table 1. SNR index for ICA and FSS optimization were calculated, for each participant and experimental phase.

To illustrate the effects of preprocessing method, time course data of the averaged ssVEP epochs recovered by FSS across trials (blue traces) are overlapped with the results of ICA (dashed green traces) in Figure 3. The outline (equivalent to the envelope) of the blue oscillations (FSS-ssVEP) overall shows more stability and less noise across time, including in the baseline segment, compared to the green (ICA) envelope, which tends to be larger in the baseline segment and shows more temporal variation during the steady-state segment. Note that data from all phases were used in the joint optimization process. This will in some cases yield a less favorable signal-to-noise ratio compared to identifying optimal components for each single phase. However, a joint optimization across phases is expected to achieve more robust and more externally valid results, being based on more trials. On the other hand, noise such as alpha oscillations (generally in the range 8–13 Hz), or the absence of an ssVEP signal, may interfere with the optimization procedure, producing divergent results as shown for participant number 07.

FIGURE 3

Figure 3. Maps of the averaged ssVEP epochs recovered by FSS (blue oscillations) along with the envelope of ssVEP optimized by ICA (dashed green line), for each participant (row) and experimental phase (column): habituation (HAB), acquisition (ACQ), and extinction (EXT). Stimulus-averaged epochs in the −1 to 7 s time period, with pre-stimulus −1 to 0 s and post-stimulus 0–7 s.

3.2. ssVEP-BOLD Fusion Results

Figure 4 and Table 2 show the ssVEP-BOLD co-variation with a standard 4 s delay. For each participant, voxels that survived individual permutation-controlled thresholding were kept and submitted to a group-level binomial test. As expected, the ssVEP-BOLD correlation map contained extended visual cortical areas, including the calcarine fissure, cuneus, occipital gyrus, and fusiform gyrus. Additional areas of ssVEP-BOLD co-variation were seen in the postcentral cortex, the rolandic operculum, and superior temporal gyrus. During the acquisition phase in which aversive learning occurs, we observed additional correlations between ssVEP and BOLD in the anterior cingulate cortex (ACC), precuneus, as well as the middle and superior temporal gyrus. Note that small but robust correlations between electrophysiological and hemodynamic measures during the same process have been consistently reported across different species, including human beings (see e.g., Boynton, 2011), for a discussion of this problem.

FIGURE 4

Figure 4. Maps of binomial test for FSS-ssVEP and BOLD fusion at each experimental phase: habituation (HAB), acquisition (ACQ), and extinction (EXT), threshold at 0.5.

TABLE 2

Table 2. Results for time-series based sVEP-BOLD correlation analysis.

To establish that using the timing information contributes to suppressing noise and increase the robustness and validity of the fusion result, we set λ to zero so that the contrast function was only driven by independency term J, allowing the ICA to use data from all three phases equally for optimization by simulated annealing. Contrasting the spatial distribution of ssVEP-fMRI correlation maps between FSS and ICA, as shown in Figure 5, that ICA yielded a less coherent, and overall noisier spatial pattern, with fewer voxels in visual areas displaying the expected effect. In addition, when comparing the scalp topographical distribution of FSS weights with ICA weights, the FSS with the focal occipital topography is physiologically more plausible than the distributed ICA topography.

FIGURE 5

Figure 5. Maps of binomial test result for FSS-ssVEP and BOLD (red) in contrast to the ICA-ssVEP and BOLD (blue), acquisition phase (ACQ). The scalp topographical distribution of FSS weights are compared with ICA.

4. Discussion

The present study aimed to illustrate the usefulness of FSS as a preprocessing step for joint analyses of fMRI-BOLD and ssVEPs, strong oscillatory signals originating in visual cortex and modulated by re-entrant inputs from a number of anterior structures (Keil et al., 2009; Norcia et al., 2015). We used data from an experiment in which participants were aversively conditioned to associate the orientation of a phase-reversing grating with a mild electric shock (Petro et al., 2017). The experiment contained three distinct phases, which were used here to examine the convergence of ssVEP-BOLD covariation maps across three re-tests, leveraging the within-participant design of this study. We found substantial overlap of ssVEP-fMRI co-variation maps in extended visual cortex, superior temporal cortex, and parietal cortex. This is consistent with the sources of ssVEP established through intracranial recordings in humans and experimental animals (see e.g., Wittevrongel et al., 2018). Furthermore, the cross-participant and re-test convergence along with the large spatial extent of the ssVEP-fMRI covariation in visual cortex observed in the present study exceeds earlier analyses of this same data-set with different methods, i.e., parametric modulation analysis as implemented in SPM (Petro et al., 2017), and multivariate correlation analyses without the FSS preprocessing step (Ji et al., 2018). Because the covariation of both signal modalities across the time series in visual cortex can be considered a gold standard for analytical strategies, the fact that ssVEP-fMRI correlations were observed reliably in extended visual cortex suggests that FSS represents a desirable preprocessing step, positively affecting the validity of EEG-fMRI fusion.

Although no explicit comparisons between experimental phases were made, the ssVEP-fMRI co-variation map during the acquisition phase, in which the visual cues were occasionally paired with electric shocks, showed co-variation of ssVEP amplitude and BOLD signal in additional parietal and anterior structures. These structures included the anterior cingulate gyrus, known to be involved in fear acquisition and associative learning in general (e.g., Sehlmeyer et al., 2009; Fullana et al., 2016). As such, this finding further supports the validity of FSS-based preprocessing for ssVEP-fMRI fusion.

It should be noted that this study was limited in terms of sample size and in terms of signal-to-noise, low in some participants, despite the precise definition of ssVEPs in the frequency domain (Norcia et al., 2015). Although information theoretic approaches (Ostwald et al., 2010, 2011) have been proposed, which apply information criteria such as mutual information to quantify the association of multimodal signals, it is difficult to estimate the probability distribution reliably with small sample size, and the results are sensitive to free parameters. Furthermore, it is not within the scope of this manuscript to quantitatively compare many alternative pipelines to establish the benefits of FSS relative to other methods. Such analyses have been performed with other data types (Porcaro et al., 2010). However, both ongoing and future research will establish and further optimize preprocessing steps for ssVEP-fMRI fusion.

5. Conclusions

In conclusion, the present study demonstrated that FSS, in conjunction with a multivariate correlation and feature selection approach proposed previously (Ji et al., 2018), yields robust ssVEP-fMRI co-variation maps with high face validity (i.e., high correspondence between BOLD and ssVEP in visual cortex, where strong responses are expected to during visual stimulation) and with high replicability across re-tests (i.e., across three experimental phases) in the same participants. The work is important to neurorobotics for several reasons: (1) This approach could be used to examine a wide range of hypotheses regarding the neurophysiology of visual cortex and its function for biological and artificial systems, (2) it provides a means to conduct the multi-modal information fusion and correlation analysis to recover the underlying neurophysiological mechanisms, and (3) the core idea of the approach are inspiring for other brain data processing and the knowledge can be transferred to other applications as well.

Ethics Statement

All procedures were approved by the institutional review board of the University of Florida and were consistent with the Declaration of Helsinki on studies with human participants.

Author Contributions

All authors had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. AK and HJ: conceptualization. HJ, ZY, and AK: methodology. HJ and NP: investigation. NP: formal analysis. HJ, AK, and BC: writing. BC and AK: supervision. BC, NZ, and AK: funding acquisition.

Funding

This work was supported by 973 Program 2015CB351703 and National Natural Science Foundation Grant 91648208 to BC and by National Institute of Mental Health grants R01 MH112558 and R01 MH097320 to AK. The funding sources had no involvement in the study design.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Acharjee, P. P., Phlypo, R., Wu, L., Calhoun, V. D., and Adalı, T. (2015). Independent vector analysis for gradient artifact removal in concurrent eeg-fmri data. IEEE Trans. Biomed. Eng. 62, 1750–1758. doi: 10.1109/TBME.2015.2403298

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, P. J., Josephs, O., and Turner, R. (2000). A method for removing imaging artifact from continuous EEG recorded during functional MRI. Neuroimage 12, 230–239. doi: 10.1006/nimg.2000.0599

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, P. J., Polizzi, G., Krakow, K., Fish, D. R., and Lemieux, L. (1998). Identification of EEG events in the MR scanner: the problem of pulse artifact and a method for its subtraction. Neuroimage 8, 229–239. doi: 10.1006/nimg.1998.0361

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, A. K. (2005). Affective influences on the attentional dynamics supporting awareness. J. Exp. Psychol. 134, 258–281. doi: 10.1037/0096-3445.134.2.258

PubMed Abstract | CrossRef Full Text | Google Scholar

Attar, C. H., Andersen, S. K., and Müller, M. M. (2010). Time course of affective bias in visual attention: convergent evidence from steady-state visual evoked potentials and behavioral data. Neuroimage 53, 1326–1333. doi: 10.1016/j.neuroimage.2010.06.074

CrossRef Full Text | Google Scholar

Barbati, G., Sigismondi, R., Zappasodi, F., Porcaro, C., Graziadio, S., and Valente, G., et al. (2006). Functional source separation from magnetoencephalographic signals. Hum. Brain Mapp. 27, 925–934. doi: 10.1002/hbm.20232

PubMed Abstract | CrossRef Full Text | Google Scholar

Bell, A. J., and Sejnowski, T. J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Comput. 7, 1129–1159. doi: 10.1162/neco.1995.7.6.1129

PubMed Abstract | CrossRef Full Text | Google Scholar

Bénar, C.-G., Schön, D., Grimault, S., Nazarian, B., Burle, B., and Roth, M., et al. (2007). Single-trial analysis of oddball event-related potentials in simultaneous EEG-FMRI. Hum. Brain Mapp. 28, 602–613. doi: 10.1002/hbm.20289

PubMed Abstract | CrossRef Full Text | Google Scholar

Boynton, G. M. (2011). Spikes, BOLD, attention, and awareness: a comparison of electrophysiological and fMRI signals in V1. J. Vis. 11:12. doi: 10.1167/11.5.12

PubMed Abstract | CrossRef Full Text | Google Scholar

Bracewell, R. N., and Bracewell, R. N. (1986). The Fourier Transform and Its Applications, Vol. 31999. New York, NY: McGraw-Hill.

Google Scholar

Buxton, R. B., Wong, E. C., and Frank, L. R. (1998). Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn. Reson. Med. 39, 855–864. doi: 10.1002/mrm.1910390602

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., He, C., Wang, Z. J., and McKeown, M. J. (2013). An IC-PLS framework for group corticomuscular coupling analysis. IEEE Trans. Biomed. Eng. 60, 2022–2033. doi: 10.1109/TBME.2013.2248059

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Peng, H., Yu, F., and Wang, K. (2017). Independent vector analysis applied to remove muscle artifacts in EEG data. IEEE Trans. Instrum. Meas. 66, 1770–1779. doi: 10.1109/TIM.2016.2608479

CrossRef Full Text | Google Scholar

Chen, X., Wang, Z. J., and McKeown, M. (2016). Joint blind source separation for neurophysiological data analysis: multiset and multimodal methods. IEEE Signal Process. Mag. 33, 86–107. doi: 10.1109/MSP.2016.2521870

CrossRef Full Text | Google Scholar

Chen, X., Xu, X., Liu, A., McKeown, M. J., and Wang, Z. J. (2018). The use of multivariate EMD and CCA for denoising muscle artifacts from few-channel EEG recordings. IEEE Trans. Instrum. Meas. 67, 359–370. doi: 10.1109/TIM.2017.2759398

CrossRef Full Text | Google Scholar

Comon, P. (1994). Independent component analysis, a new concept? Signal Process. 36, 287–314. doi: 10.1016/0165-1684(94)90029-9

CrossRef Full Text | Google Scholar

Debener, S., Ullsperger, M., Siegel, M., and Engel, A. K. (2006). Single-trial EEG–fMRI reveals the dynamics of cognitive function. Trends Cogn. Sci. 10, 558–563. doi: 10.1016/j.tics.2006.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Debener, S., Ullsperger, M., Siegel, M., Fiehler, K., Von Cramon, D. Y., and Engel, A. K. (2005). Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identifies the dynamics of performance monitoring. J. Neurosci. 25, 11730–11737. doi: 10.1523/JNEUROSCI.3286-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Russo, F., Pitzalis, S., Aprile, T., Spitoni, G., Patria, F., and Stella, A., et al. (2007). Spatiotemporal analysis of the cortical sources of the steady-state visual evoked potential. Hum. Brain Mapp. 28, 323–334. doi: 10.1002/hbm.20276

PubMed Abstract | CrossRef Full Text | Google Scholar

Eichele, T., Specht, K., Moosmann, M., Jongsma, M. L., Quiroga, R. Q., and Nordby, H., et al. (2005). Assessing the spatiotemporal evolution of neuronal activation with single-trial event-related potentials and functional MRI. Proc. Natl. Acad. Sci. U.S.A. 102, 17798–17803. doi: 10.1073/pnas.0505508102

PubMed Abstract | CrossRef Full Text | Google Scholar

Eickhoff, S. B., Stephan, K. E., Mohlberg, H., Grefkes, C., Fink, G. R., and Amunts, K., et al. (2005). A new spm toolbox for combining probabilistic cytoarchitectonic maps and functional imaging data. Neuroimage 25, 1325–1335. doi: 10.1016/j.neuroimage.2004.12.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Falotico, E., Vannucci, L., Ambrosano, A., Albanese, U., Ulbrich, S., and Vasquez Tieck, J. C., et al. (2017). Connecting artificial brains to robots in a comprehensive simulation framework: the neurorobotics platform. Front. Neurorobot. 11:2. doi: 10.3389/fnbot.2017.00002

PubMed Abstract | CrossRef Full Text | Google Scholar

Fullana, M., Harrison, B., Soriano-Mas, C., Vervliet, B., Cardoner, N., and Àvila-Parcet, A., et al. (2016). Neural signatures of human fear conditioning: an updated and extended meta-analysis of fMRI studies. Mol. Psychiatry 21:500. doi: 10.1038/mp.2015.88

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrmann, C. S., and Debener, S. (2008). Simultaneous recording of EEG and BOLD responses: a historical perspective. Int. J. Psychophysiol. 67, 161–168. doi: 10.1016/j.ijpsycho.2007.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, N. E. (2014). Hilbert-Huang Transform and Its Applications, Vol. 16. Singapore: World Scientific.

Google Scholar

Hyvarinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Netw. 10, 626–634. doi: 10.1109/72.761722

PubMed Abstract | CrossRef Full Text | Google Scholar

Hyvärinen, A., and Oja, E. (2000). Independent component analysis: algorithms and applications. Neural Netw. 13, 411–430. doi: 10.1016/S0893-6080(00)00026-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ingber, L. (2000). Adaptive simulated annealing (ASA): lessons learned. arXiv [Preprint]. arXiv:cs/0001018.

Google Scholar

James, C. J., and Hesse, C. W. (2004). Independent component analysis for biomedical signals. Physiol. Meas. 26:R15. doi: 10.1007/978-3-540-30110-3_129

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, H., Petro, N. M., Chen, B., Yuan, Z., Wang, J., and Zheng, N., et al. (2018). Cross multivariate correlation coefficients as screening tool for analysis of concurrent EEG-fMRI recordings. J. Neurosci. Res. 96, 1159–1175. doi: 10.1002/jnr.24217

PubMed Abstract | CrossRef Full Text | Google Scholar

Keil, A., Gruber, T., Müller, M. M., Moratti, S., Stolarova, M., and Bradley, M. M., et al. (2003). Early modulation of visual perception by emotional arousal: evidence from steady-state visual evoked brain potentials. Cogn. Affect. Behav. Neurosci. 3, 195–206. doi: 10.3758/CABN.3.3.195

PubMed Abstract | CrossRef Full Text | Google Scholar

Keil, A., Miskovic, V., Gray, M. J., and Martinovic, J. (2013). Luminance, but not chromatic visual pathways, mediate amplification of conditioned danger signals in human visual cortex. Eur. J. Neurosci. 38, 3356–3362. doi: 10.1111/ejn.12316

CrossRef Full Text | Google Scholar

Keil, A., Sabatinelli, D., Ding, M., Lang, P. J., Ihssen, N., and Heim, S. (2009). Re-entrant projections modulate visual cortex in affective perception: evidence from granger causality analysis. Hum. Brain Mapp. 30, 532–540. doi: 10.1002/hbm.20521

PubMed Abstract | CrossRef Full Text | Google Scholar

Keil, A., Smith, J. C., Wangelin, B. C., Sabatinelli, D., Bradley, M. M., and Lang, P. J. (2008). Electrocortical and electrodermal responses covary as a function of emotional arousal: a single-trial analysis. Psychophysiology 45, 516–523. doi: 10.1111/j.1469-8986.2008.00667.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P. (1983). Optimization by simulated annealing. Science 220, 671–680. doi: 10.1126/science.220.4598.671

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Huang, H., McGinnis-Deweese, M., Keil, A., and Ding, M. (2012). Neural substrate of the late positive potential in emotional processing. J. Neurosci. 32, 14563–14572. doi: 10.1523/JNEUROSCI.3109-12.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Logothetis, N. K. (2015). Neural-event-triggered fMRI of large-scale neural networks. Curr. Opin. Neurobiol. 31, 214–222. doi: 10.1016/j.conb.2014.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Logothetis, N. K., Pauls, J., Augath, M., Trinath, T., and Oeltermann, A. (2001). Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 150–157. doi: 10.1038/35084005

PubMed Abstract | CrossRef Full Text | Google Scholar

Makeig, S., Westerfield, M., Jung, T.-P., Enghoff, S., Townsend, J., and Courchesne, E., et al. (2002). Dynamic brain sources of visual evoked responses. Science 295, 690–694. doi: 10.1126/science.1066168

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantini, D., Perrucci, M. G., Cugini, S., Ferretti, A., Romani, G. L., and Del Gratta, C. (2007). Complete artifact removal for EEG recorded during continuous fMRI using independent component analysis. Neuroimage 34, 598–607. doi: 10.1016/j.neuroimage.2006.09.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Mulert, C., Seifert, C., Leicht, G., Kirsch, V., Ertl, M., and Karch, S., et al. (2008). Single-trial coupling of EEG and fMRI reveals the involvement of early anterior cingulate cortex activation in effortful decision making. Neuroimage 42, 158–168. doi: 10.1016/j.neuroimage.2008.04.236

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, M. M., Teder, W., and Hillyard, S. A. (1997). Magnetoencephalographic recording of steady-state visual evoked cortical activity. Brain Topogr. 9, 163–168. doi: 10.1007/BF01190385

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, M. M., Teder-Salejarvi, W., and Hillyard, S. A. (1998). The time course of cortical facilitation during cued shifts of spatial attention. Nat. Neurosci. 1, 631–634. doi: 10.1038/2865

PubMed Abstract | CrossRef Full Text | Google Scholar

Naik, G. R., and Wang, W. (2014). Blind Source Separation: Advances in Theory Algorithms and Applications. Signals and Communication Technology Series. Berlin: Springer.

Google Scholar

Norcia, A. M., Appelbaum, L. G., Ales, J. M., Cottereau, B. R., and Rossion, B. (2015). The steady-state visual evoked potential in vision research: a review. J. Vis. 15:4. doi: 10.1167/15.6.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez, P. L., and Srinivasan, R. (2006). Electric Fields of the Brain: The Neurophysics of EEG. New York, NY: Oxford University Press.

Google Scholar

Oess, T., Krichmar, J. L., and Röhrbein, F. (2017). A computational model for spatial navigation based on reference frames in the hippocampus, retrosplenial cortex, and posterior parietal cortex. Front. Neurorobot. 11:4. doi: 10.3389/fnbot.2017.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Öhman, A., and Soares, J. J. (1998). Emotional conditioning to masked stimuli: expectancies for aversive outcomes following nonrecognized fear-relevant stimuli. J. Exp. Psychol. 127:69. doi: 10.1037//0096-3445.127.1.69

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostwald, D., Porcaro, C., and Bagshaw, A. P. (2010). An information theoretic approach to EEG–fMRI integration of visually evoked responses. Neuroimage 49, 498–516. doi: 10.1016/j.neuroimage.2009.07.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostwald, D., Porcaro, C., and Bagshaw, A. P. (2011). Voxel-wise information theoretic EEG-fMRI feature integration. Neuroimage 55, 1270–1286. doi: 10.1016/j.neuroimage.2010.12.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Penny, W. D., Friston, K. J., Ashburner, J. T., Kiebel, S. J., and Nichols, T. E. (2011). Statistical Parametric Mapping: The Analysis of Functional Brain Images. London: Elsevier.

Google Scholar

Petro, N. M., Gruss, L. F., Yin, S., Huang, H., Miskovic, V., and Ding, M., et al. (2017). Multimodal imaging evidence for a frontoparietal modulation of visual cortex during the selective processing of conditioned threat. J. Cogn. Neurosci. 29, 953–967. doi: 10.1162/jocn_a_01114

PubMed Abstract | CrossRef Full Text | Google Scholar

Porcaro, C., Ostwald, D., and Bagshaw, A. P. (2010). Functional source separation improves the quality of single trial visual evoked potentials recorded during concurrent EEG-fMRI. Neuroimage 50, 112–123. doi: 10.1016/j.neuroimage.2009.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Porcaro, C., Ostwald, D., Hadjipapas, A., Barnes, G. R., and Bagshaw, A. P. (2011). The relationship between the visual evoked potential and the gamma band investigated by blind and semi-blind methods. NeuroImage 56, 1059–1071. doi: 10.1016/j.neuroimage.2011.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Regan, D. (1989). Human Brain Electrophysiology: Evoked Potentials and Evoked Magnetic Fields in Science and Medicine. New York, NY: Elsevier.

Google Scholar

Sammer, G., Blecker, C., Gebhardt, H., Kirsch, P., Stark, R., and Vaitl, D. (2005). Acquisition of typical EEG waveforms during fMRI: SSVEP, LRP, and frontal theta. Neuroimage 24, 1012–1024. doi: 10.1016/j.neuroimage.2004.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Sehlmeyer, C., Schöning, S., Zwitserlood, P., Pfleiderer, B., Kircher, T., and Arolt, V., et al. (2009). Human fear conditioning and extinction in neuroimaging: a systematic review. PLoS ONE 4:e5865. doi: 10.1371/journal.pone.0005865

PubMed Abstract | CrossRef Full Text | Google Scholar

Serences, J. T. (2004). A comparison of methods for characterizing the event-related bold timeseries in rapid fmri. Neuroimage 21, 1690–1700. doi: 10.1016/j.neuroimage.2003.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Sokolov, E. N. (1963). Perception and the Conditioned Reflex. New York, NY: The Macmillan Company.

Google Scholar

Spekreijse, H., Dagnelie, G., Maier, J., and Regan, D. (1985). Flicker and movement constituents of the pattern reversal response. Vis. Res. 25, 1297–1304. doi: 10.1016/0042-6989(85)90045-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Todd, R., and Thompson, E. (2015). Strengthening emotion-cognition integration. Behav. Brain Sci. 38:e87. doi: 10.1017/S0140525X14001022

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C.-C., Thorpe, C., Thrun, S., Hebert, M., and Durrant-Whyte, H. (2007). Simultaneous localization, mapping and moving object tracking. Int. J. Robot. Res. 26, 889–916. doi: 10.1177/0278364907081229

CrossRef Full Text | Google Scholar

Wang, J., and Zheng, N. (2014). Measures of linear correlation for multiple variables. arXiv [Preprint]. arXiv:1401.4827.

Google Scholar

Wieser, M. J., and Keil, A. (2011). Temporal trade-off effects in sustained attention: dynamics in visual cortex predict the target detection performance during distraction. J. Neurosci. 31, 7784–7790. doi: 10.1523/JNEUROSCI.5632-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wittevrongel, B., Khachatryan, E., Hnazaee, M. F., Carrette, E., De Taeye, L., and Meurs, A., et al. (2018). Representation of steady-state visual evoked potentials elicited by luminance flicker in human occipital cortex: an electrocorticography study. Neuroimage 175, 315–326. doi: 10.1016/j.neuroimage.2018.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar