Simultaneous EEG-fMRI: Denoising, Frequency Analysis, and Correlation Mapping for Pediatric Clinic Planning

Assignment 3: Simultaneous EEG-fMRI Due Friday July 7th at 11:59 PM In this assignment you will get familiar with digital signal processing including common denoising techniques (ICA), filtering, resampling, and convolution. The goal of this assignment is to denoise an EEG dataset acquired during a simultaneous EEG-fMRI experiment, extract some frequency-domain measures from the denoised EEG signal, and correlate with fMRI image (which is a 4d image, each voxel is a time series). The end-result will be whole-brain correlation maps showing which brain areas correlate best with EEG. Background: simultaneous EEG-fMRI is a multimodal neuroimaging technique that combines EEG with fMRI, yielding two separate datasets, which when combined achieve excellent temporal resolution (due to the high temporal sampling rate of EEG) and excellent spatial resolution (due to the high spatial resolution of fMRI). Simultaneous EEG-fMRI allows to investigate the precise brain areas involved in the generation of frequency-specific brain oscillations. Here, we are interested in which brain areas correlate with alpha (8-13 Hz) and gamma (40-70 Hz) band power. Therefore, we need to get the power from the EEG data and correlate it with BOLD signal in each voxel. To do so, there are some basic preprocessing steps we need to perform: Step 1: Understand the EEG and fMRI datasets (LINK). The dataset contains 3 subjects and 3 scans per subject (rest, gamma_01, and gamma_02), each scan has an EEG and an fMRI dataset. This dataset features a resting state scan (eyes closed, relaxed) and two visual experiments (watching a screen with some visual stimuli). The EEG dataset contains 64 separate time series (one time series per electrode). The 32nd electrode (index 31) is an ECG (electrocardiogram) electrode which records from the heart, the other 63 electrodes record from the brain. The sampling rate of the EEG is 5000 Hz. The fMRI dataset is a 4D image (3D + time). Information about spatial and temporal resolution of the fMRI images is contained in the .nii headers. Anatomical images (T1 and SWI) are also included. Step 2: denoise the EEG dataset. The EEG data was recorded inside the MRI scanner (in order to be acquired simultaneously with fMRI), so it has two major artifacts 1) the gradient artifact (figure 1) and 2) the ballistocardiogram artifact (figure 2). The gradient artifact must be removed first, the ballistocardiogram artifact can be removed second.

a) Load EEG dataset into python using mne b) Remove gradient and ballistocardiogram artifacts (figure 3)

a. Use average artifact subtraction for both (see pdf on moodle for description of this) c) Run ICA (figure 5) on EEG signal to isolate good and bad components

a. can use the mne implementation of fastICA d) Find the ‘good’ and ‘bad’ ICA components by visualizing the topography and power spectrum

a. Identifying ‘good’ and ‘bad’ components can be tricky and requires some experience. I will discuss this in future lectures.

e) Bandpass filter the denoised EEG data to separate it into two separate frequency bands: 1. Alpha band (8-13 Hz) 2. Gamma band (40-70 Hz) Once step (e) is done, rectify (take absolute value) the bandpass filtered signals and average across all ‘good’ components.

Step 3: denoise the fMRI dataset. the fMRI dataset must be motion corrected and bandpass filtered:

a) Load fMRI dataset into python using nibabel b) Motion correct by registering every volume to the first volume

a. find some python library that can do this, or use FSL, AFNI or ANTs c) Apply 0.01 – 0.1 Hz bandpass filter to the time series in each voxel

a. This will remove all frequencies outside 0.01 – 0.1 Hz. The remaining signal should be ‘good’ in the sense that it will not contain much noise (most fMRI noise is outside the 0.01-0.1 Hz. band)

Step 4: combine EEG and fMRI datasets. After completing the pre-processing for both datasets (step 2,3) it is time to combine the datasets and get our correlation maps:

a) Resample the smooth, rectified bandpass filtered EEG component to match the fMRI temporal resolution.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3351136/
https://drive.google.com/file/d/1HsjKiG32lOfmwcAdNOgD2AN0LjJzI8DN/view?usp=drive_link
https://mne.tools/stable/index.html
b) Convolve the resampled, bandpass-filtered EEG power with the canonical hemodynamic response function. This will give the ‘expected’ fMRI response to the frequency-specific brain activity.

c) Correlate the convolved EEG alpha (8-13 Hz) and gamma (40-70 Hz) power with fMRI in each voxel, yielding two separate 3d images where the value in each voxel is a correlation coefficient.

i. (optional) It might also help tp bandpass filter the rectified, resampled EEG alpha between 0.01 – 0.1 Hz, because power can drift over time, inducing low-frequency power that can distort the correlation.

Step 5: combine subjects and scans for final result. At the end, you will have 6 correlation maps per subject (3 scans x 2 frequency bands). The final result is the grand average correlation map (across the 3 subjects) for each frequency band. Show the following result in your pdf report: Plot 1-3: show the correlation map for each subject (two maps per subject) overlaid on each subject’s T1-weighted image. Average across all 3 scans (resting state + gamma x2) when showing these plots. Plot 4: show the same as plot 1-3, but average the correlation maps across all subjects (use russell’s T1 as the underlay). You will need to register Alex and Genevieve’s T1 images into Russell space and then apply the matrix to the correlation map, so all 3 subjects are in the same space before averaging.

You can show a single axial slice in each plot. Choose the slice that shows the visual cortex the best (it’s usually about 10- 15 mm above the cerebellum). Other results to show in your pdf: besides the final result from step 5, you should show the ICA components from a single subject, the EEG signal from a single channel after a) gradient artifact removal b) BCG artifact removal and c) ‘bad’ ICA component subtraction. I will release video lectures covering all the steps in the next week. Figure 1: Before gradient artifact subtraction. The raw EEG signal is heavily corrupted by the MRI scanner gradient pulse.

Figure 2: After gradient artifact subtraction but before ballistocardiogram artifact subtraction. Pulse artifact visible.

Figure 3: comparing before (blue) and after (orange) ballistocardiogram artifact subtraction

Figure 4: an example of detected ballistocardiogram artifact peaks (scipy.signal.find_peaks).

Figure 5: ICA component weight maps for components 0-19 (out of a total 63, only 0-19 shown here). Components 0,1 and 2 are definitely artifacts (Eye blinks and residual bcg). Components 5,15, and 19 are definitely neural in origin (what we want). The other components its more difficult to tell, usually it’s best to discard components when unsure.