Inter-Subject Correlation and Inter-Subject Functional Correlation

Contributions

The functional connectivity methods that we used in previous notebooks compared time series of BOLD activity between voxels within participant to infer how different regions of the brain were interacting. However, BOLD activity contains multiple components (Figure a):

  1. Task-based/stimulus-evoked signal that is reliable across participants
  2. Intrinsic fluctuations in neural activity that are participant specific
  3. Scanner or physiological noise

In this notebook, we will consider methods that combine data across participants to eliminate #2 and #3 when calculating fMRI reliability (intersubject correlation, ISC, Hasson et al., 2004) and connectivity (intersubject functional correlation, ISFC, Simony et al., 2016). ISC and ISFC help isolate #1 because it is the only component that ought to be shared across participants.

Figure b,c show how ISC differs from functional connectivity: rather than correlating brain regions, which preserves participant-specific activity and noise, ISC correlates between the brains of different participants in order to capture only the activity that is shared. In ISC, this correlation is done for every voxel in the brain to the matching voxel in other brains, producing a full brain map. Figure e shows this as the diagonal of a correlation matrix, where each cell corresponds to a voxel in subject X correlated with the same anatomical voxel in subject Y. In practice, to simplify the computation and the interpretation it is typical for ISC to compare each individual participant with the average of all other participants.

Figure d shows ISFC: the correlation of every voxel in one participant with every other voxel in another participant (or average of other participants). This is like FCMA except it is between participants rather than within participants. In fact, these analyses use the same computational tricks. ISFC is valuable because it allows us to identify activity coupling in voxels that are not aligned across participants: the off diagonal in Figure e represents correlations for voxels in different parts of the brain.

alt text

We will use ISC and ISFC to identify brain regions that respond preferentially to narrative stories, rather than to a random assortment of words (replicating Simony et al., 2016). Furthermore, seed-based connectivity analysis does not show differences between resting state, random words, and intact narratives, but ISFC does distinguish between these conditions (Simony et al., 2016). Thus, ISFC shows greater sensitivity to the task than seed-based functional connectivity.

Goal of this script

1. To run intersubject correlation (ISC).
2. To run intersubject functional correlation (ISFC).  
3. Use ISFC to examine how a network of brain regions that respond to narrative stimuli.  


Table of Contents

1. The ISC-ISFC Workflow

2. ISC

2.1 The "Pieman" data
2.2 Data file preparation
2.3 Compute ISC
2.4 ISC with statistical tests

3. ISFC

3.1 Parcel the data
3.2 Compute FC and ISFC

4. Spatial Correlation

4.1 Spatial inter-subject correlation

Exercises

1 2 3 4 5 6 7 8 9
Novel contribution

1. The ISC-ISFC workflow

The following sequence of steps are recommended for successfully running ISC and ISFC using BrainIAK.

  1. Data Preparation: Organize a data directory with fMRI subject data that you want to process. All subjects must be in the same anatomical space for analysis. Also you need to create a whole-brain mask. The outcome of this is an array of anatomically-aligned and temporally-aligned brain data.

  2. Compute ISC: The ISC function computes correlations across subjects for corresponding voxels in the mask. It uses the compute_correlation function in BrainIAK, which is optimized for fast execution (and was used in FCMA).

  3. Permutation Test for ISC: Perform statistical analysis to determine significant correlation values for ISC.

  4. Compute ISFC: The ISFC function computes correlations for every voxel in one subject with every other voxel averaged across subjects.

  5. Cluster the ISFC results: Create clusters based on the correlation values.

  6. Perform ISFC permutation: Perform permutation tests to determine the significance of the results.

In [1]:
import warnings
import sys 
if not sys.warnoptions:
    warnings.simplefilter("ignore")
import os 
import glob
import time
from copy import deepcopy
import numpy as np
import pandas as pd 

from nilearn import datasets
from nilearn import surface
from nilearn import plotting
from nilearn.input_data import NiftiMasker, NiftiLabelsMasker
import nibabel as nib

from brainiak import image, io
from brainiak.isc import isc, isfc, permutation_isc
import matplotlib.pyplot as plt
import seaborn as sns 

%autosave 5
%matplotlib inline
sns.set(style = 'white', context='talk', font_scale=1, rc={"lines.linewidth": 2})
Autosaving every 5 seconds

2. ISC

2.1 The "Pieman" data

For this script we will use the "Pieman" dataset from Simony et al. (2016). A description of the dataset is as follows:

18 native English speakers were scanned (15 females, ages: 18–31), corresponding to the replication dataset from the Pieman study.
Stimuli for the experiment were generated from a 7 min real life story ("Pie Man", Jim O'Grady) recorded at a live storytelling performance ("The Moth" storytelling event, New York City). Subjects listened to the story from beginning to end (intact condition). In addition, subjects listened to scrambled versions of the story, which were generated by dividing the original stimulus into segments of different timescales (paragraphs and words) and then permuting the order of these segments. To generate the scrambled stimuli, the story was segmented manually by identifying the end points of each word and paragraph. Two adjacent short words were assigned to a single segment in cases where we could not separate them. Following segmentation, the intact story was scrambled at two timescales: short—‘words’ (W; 608 words, 0.7±0.5 s each) and long—‘paragraphs’ (P; 11 paragraphs, 38.1±17.6 s each). Laughter and applause were classified as single word events (4.4% of the words). Twelve seconds of neutral music and 3 s of silence preceded, and 15 s of silence followed, each playback in all conditions. These music and silence periods were discarded from all analyses.

More details about the experiment may be accessed in the methods section of the paper.

2.2 Data File Preparation

Loading and preparing the data:

BrainIAK has methods to efficiently load data. We have used some of these functions in previous notebooks.

load_images: reads data from all subjects in a list that you provide. This is like the function load_images_from_dir but here we specify the names manually.

load_boolean_mask: Create a binary mask from a brain volume

mask_images: Loads the brain images and masks them with the mask provided

image.MaskedMultiSubjectData.from_masked_images: Creates a list of arrays, with each item in the list corresponding to one subject's data. This data format is accepted by the BrainIAK ISC and ISFC function.

In [2]:
# Set up experiment metadata

from utils import pieman2_dir, results_path
print('Data directory is: %s' % pieman2_dir)

dir_mask = os.path.join(pieman2_dir, 'masks/')
mask_name = os.path.join(dir_mask, 'avg152T1_gray_3mm.nii.gz')
all_task_names = ['word', 'intact1']
all_task_des = ['word level scramble', 'intact story']
n_subjs_total = 18
group_assignment_dict = {task_name: i for i, task_name in enumerate(all_task_names)}

# Where do you want to store the data
dir_out = results_path + 'isc/'
if not os.path.exists(dir_out):
    os.makedirs(dir_out)
    print('Dir %s created ' % dir_out)
Data directory is: /jukebox/norman/mkumar/tutorial_data/brainiak_datasets/Pieman2

Helper functions

We provide helper functions to load the data.

Memory limits Be aware this is going to be run on 18 participants and may push the limits of your memory and computational resources if you are on a laptop. If you want to run it on fewer participants to protect memory, change `n_subjs` to be lower (e.g. 10); however, the anticipated results may not generalize to lower sample sizes.
In [3]:
# Reduce the number of subjects per condition to make this notebook faster 
upper_limit_n_subjs = 18

def get_file_names(data_dir_, task_name_, verbose = False):
    """
    Get all the participant file names
    
    Parameters
    ----------
    data_dir_ [str]: the data root dir
    task_name_ [str]: the name of the task 
    
    Return
    ----------
    fnames_ [list]: file names for all subjs
    """
    c_ = 0 
    fnames_ = []
    # Collect all file names 
    for subj in range(1, n_subjs_total): 
        fname = os.path.join(
            data_dir_, 'sub-%.3d/func/sub-%.3d-task-%s.nii.gz' % (subj, subj, task_name_))
        
        # If the file exists
        if os.path.exists(fname):
            
            # Add to the list of file names 
            fnames_.append(fname)
            if verbose: 
                print(fname)
            c_+= 1
            if c_ >= upper_limit_n_subjs: 
                break
    return fnames_
In [4]:
"""load brain template"""

# Load the brain mask
brain_mask = io.load_boolean_mask(mask_name)

# Get the list of nonzero voxel coordinates
coords = np.where(brain_mask)

# Load the brain nii image
brain_nii = nib.load(mask_name)
In [5]:
"""load bold data"""

# load the functional data 
fnames = {}
images = {}
masked_images = {}
bold = {}
group_assignment = []
n_subjs = {}

for task_name in all_task_names: 
    fnames[task_name] = get_file_names(pieman2_dir, task_name)
    images[task_name] = io.load_images(fnames[task_name]) 
    masked_images[task_name] = image.mask_images(images[task_name], brain_mask) 
    # Concatenate all of the masked images across participants  
    bold[task_name] = image.MaskedMultiSubjectData.from_masked_images(
        masked_images[task_name], len(fnames[task_name])
    )
    # Convert nans into zeros
    bold[task_name][np.isnan(bold[task_name])] = 0
    # compute the group assignment label 
    n_subjs_this_task = np.shape(bold[task_name])[-1]
    group_assignment += list(
        np.repeat(group_assignment_dict[task_name], n_subjs_this_task)
    )
    n_subjs[task_name] = np.shape(bold[task_name])[-1]
    print('Data loaded: {} \t shape: {}' .format(task_name, np.shape(bold[task_name])))
Data loaded: word 	 shape: (300, 98508, 9)
Data loaded: intact1 	 shape: (300, 98508, 17)

Exercise 1: Inspect the data and report on the following details.

  • Brain template
    • Report the shape of brain_nii, brain_mask
    • Visualize brain_nii and brain_mask by plotting the 30th slice along the Z dimension.
    • Describe what coords refers to
    • Visualize coords with a 3d plot. For this, only plot every 10th point, otherwise the plot will be slow to load.
  • Brain data
    • Inspect the shape of bold. How many subjects do we have for each task condition? Do different subjects have the same number of TRs/voxels?
In [6]:
# Insert code below

2.3 Compute ISC

ISC is the correlation of each voxel's time series for a participant with the corresponding (anatomically aligned) voxel time series in the average of the other participants' brains. BrainIAK has functions for computing ISC by feeding in the concatenated participant data.

This will take about 10 minutes to complete.

In [7]:
# run ISC, loop over conditions 
isc_maps = {}
for task_name in all_task_names:
    isc_maps[task_name] = isc(bold[task_name], pairwise=False)
    print('Shape of %s condition:' % task_name, np.shape(isc_maps[task_name]))
Shape of word condition: (9, 98508)
Shape of intact1 condition: (17, 98508)

The output of ISC is a voxel by participant matrix (showing the result of each individual with the group). Below we will visualize the ISC matrix for one participant and condition back on to the brain to see where activity is correlated between participants.

In [8]:
# set params 
subj_id = 0
task_name = 'intact1'
save_data = False

# Make the ISC output a volume
isc_vol = np.zeros(brain_nii.shape)
# Map the ISC data for the first participant into brain space
isc_vol[coords] = isc_maps[task_name][subj_id, :]
# make a nii image of the isc map 
isc_nifti = nib.Nifti1Image(isc_vol, brain_nii.affine, brain_nii.header)

# Save the ISC data as a volume
if save_data: 
    isc_map_path = os.path.join(dir_out, 'ISC_%s_sub%.2d.nii.gz' % (task_name, subj_id))
    nib.save(isc_nifti, isc_map_path)
In [9]:
# Plot the data as a statmap
threshold = .2

f, ax = plt.subplots(1,1, figsize = (12, 5))
plotting.plot_stat_map(
    isc_nifti, 
    threshold=threshold, 
    axes=ax
)
ax.set_title('ISC map for subject {}, task = {}' .format(subj_id,task_name)) 
Out[9]:
Text(0.5, 1.0, 'ISC map for subject 0, task = intact1')

Exercise 2: Visualize the averaged ISC map (averaged across participants) for each task.

  • Make the averaged ISC map for the two conditions (intact, word-scrambled).
  • Visualize them using plotting.plot_stat_map.

Make sure to compare the two maps using the same xyz cut, threshold and vmax.

In [10]:
# Insert code here

This analysis was performed in volumetric space; however, nilearn makes it easy to compare this data in surface space (assuming the alignment to MNI standard is excellent). Here's an example of surface plot.

In [11]:
# set some plotting params 
subj_id = 0 
task_name = 'intact1'
threshold = .2 
view = 'medial'

# get a surface
fsaverage = datasets.fetch_surf_fsaverage5()

# Make the ISC output a volume
isc_vol = np.zeros(brain_nii.shape)
# Map the ISC data for the first participant into brain space
isc_vol[coords] = isc_maps[task_name][subj_id, :]
# make a nii image of the isc map 
isc_intact_1subj = nib.Nifti1Image(isc_vol, brain_nii.affine, brain_nii.header)

# make "texture" 
texture = surface.vol_to_surf(isc_intact_1subj, fsaverage.pial_left) 

# plot 
title_text = ('Avg ISC map, {} for one participant'.format(task_name))
surf_map = plotting.plot_surf_stat_map(
    fsaverage.infl_left, texture, 
    hemi='left', view=view, 
    title= title_text, 
    threshold=threshold, cmap='RdYlBu_r', 
    colorbar=True,
    bg_map=fsaverage.sulc_left)
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/ipykernel_launcher.py:8: VisibleDeprecationWarning: fetch_surf_fsaverage5 has been deprecated and will be removed in a future release. Use fetch_surf_fsaverage(mesh='fsaverage5')
  

Exercise 3: Visualize the averaged ISC map using surface plot.

  • Visualize the average ISC maps using plotting.plot_surf_stat_map for:
    • both conditions
    • both medial view and lateral views

Make sure you are using the same threshold and vmax for all plots.

In [12]:
# Insert code here

Exercise 4: Compare the averaged ISC map for the two task conditions. What are some brain regions showing stronger correlation in the intact story condition (vs. the word-level scramble condition)? What does this tell us about the processing of language?

Hint: The following paper this work comes from will help.

A:

2.4 ISC with statistical tests

BrainIAK provides several nonparametric statistical tests for ISC analysis (Nastase et al., 2019). Nonparametric tests are preferred due to the inherent correlation structure across ISC values—each subject contributes to the ISC of other subjects, violating assumptions of independence required for standard parametric tests (e.g., t-test, ANOVA). We will use the permutation test below.

2.4.1 Permutation test

Permutation tests are used to compute a null distribution of values. We have used permutation tests in previous notebooks and the steps outlined here are similar to what was done in prior notebooks with one small change, incorporating the group of subjects to compute the ISC:

  1. Prepare the data. Here we have two conditions (intact and word_scramble), so we compute ISC for both conditions and concatenate the data for these two conditions for all subjects.

    We use leave-one-subject-out (pairwise=False) to compute ISC and then use these correlations isc_maps_all_tasks to compute statistics.

  2. We are going to permute the condition label for each subject to simulate the randomization of conditions. To do this, we first need to assign subjects to the correct experimental conditions that they were in. We have prepared such a list of assignments when we loaded the data and stored the information in the variable: group_assignment.

  3. The next steps are executed internally in BrainIAK in the function permutation_isc:

    • For each permutation iteration:
      • BrainIAK permutes the group assignment for each subject.
      • A mean of the ISC values is then computed for this shuffled group for each condition.
      • A difference of group means is computed between each condition.
    • The difference values for all iterations is collected and forms the null distribution.
  4. Finally, we pick a threshold value that corresponds to p percent of this distribution.

permutation_isc returns the actual observed ISC values, p-values, and optionally the resampling distribution.

In [13]:
# Concatenate ISCs from both tasks
isc_maps_all_tasks = np.vstack([isc_maps[task_name] for
                                task_name in all_task_names])

print('group_assignment: {}'.format(group_assignment))
print('isc_maps_all_tasks: {}' .format(np.shape(isc_maps_all_tasks)))
group_assignment: [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
isc_maps_all_tasks: (26, 98508)
In [14]:
# permutation testing
n_permutations = 1000
summary_statistic='mean'

observed, p, distribution = permutation_isc(
    isc_maps_all_tasks, 
    pairwise=False,
    group_assignment=group_assignment, 
    summary_statistic=summary_statistic,
    n_permutations=n_permutations
)

p = p.ravel()
observed = observed.ravel()

print('observed:{}'.format(np.shape(observed)))
print('p:{}'.format(np.shape(p)))
print('distribution: {}'.format(np.shape(distribution)))
observed:(98508,)
p:(98508,)
distribution: (1000, 98508)

Exercise 5: Interpret the results from the permutation test.

  • What's the logic of the permutation test?
  • What are the outputs observed, p, distribution
  • Visualize the correlation contrast map (e.g. intact > word-level scramble) with a significance criterion (e.g. p < .005). Which region(s) showed higher ISC under the chosen contrast?

A:

3. ISFC

The goal of ISFC is to find coupling between brain regions across participants. For example the angular gyrus in subject 1 could be correlated to the pre-frontal cortex in subject 2, if they share some cognitive state. For completely random cognitive states across these two subjects, the correlation should be zero. ISFC helps us identify such commonalities across subjects.

In this section, we will compare functional connectivity vs. ISFC on the Pieman data. Whereas FC is computed within individuals, ISFC is computed between individuals. Hence the only correlations that should be robust in ISFC are those that are present across individuals. At the end of the exercises, you will qualitatively replicate Simony et al. (2016), showing that ISFC is sensitive to the cognitive state of the participants.

3.1 Parcel the data

ISFC in voxel space is very computationally intensive, so for this notebook we will divide the brain into a smaller number of parcels. We are going to use predefined ROI masks to select the voxels.

In [15]:
# load a parcel
atlas = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-2mm', symmetric_split=True)
plotting.plot_roi(atlas.maps, title='the harvard-oxford parcel')

n_regions = len(atlas.labels)-1 # rm background region 
n_TRs = np.shape(bold[task_name])[0]

print('number of voxels:\t {}'.format(np.shape(bold[task_name][1])))
print('number of parcels:\t {}'.format(n_regions))
number of voxels:	 (98508, 17)
number of parcels:	 96

Convert the bold data into ROI parcels

In [16]:
# Get a masker for the atlas 
masker_ho = NiftiLabelsMasker(labels_img=atlas.maps)

# Transform the data to the parcel space
bold_ho = {
    task_name:np.zeros((n_TRs, n_regions, n_subjs[task_name])) 
    for task_name in all_task_names}

# Collect all data 
row_has_nan = np.zeros(shape=(n_regions,), dtype=bool)
for task_name in all_task_names:
    for subj_id in range(n_subjs[task_name]):
        
        # get the data for task t, subject s 
        nii_t_s = nib.load(fnames[task_name][subj_id])
        bold_ho[task_name][:,:,subj_id] = masker_ho.fit_transform(nii_t_s)
        
        # figure out missing rois
        row_has_nan_ = np.any(np.isnan(bold_ho[task_name][:,:,subj_id]),axis=0)
        row_has_nan[row_has_nan_] = True                

# Figure out which ROI has missing values
roi_select = np.logical_not(row_has_nan)
n_roi_select = np.sum(roi_select)
rois_filtered = np.array(atlas.labels[1:])[roi_select]
bold_ho_filtered = {
    task_name:np.zeros((n_TRs, n_roi_select, n_subjs[task_name])) 
    for task_name in all_task_names
}

# Remove ROIs with missing values
for task_name in all_task_names:
    for subj_id in range(n_subjs[task_name]):
        bold_ho_filtered[task_name][:,:,subj_id] = bold_ho[task_name][:,roi_select,subj_id]
              
print('ROI selected\n {}'.format(rois_filtered))
print('ROI removed due to missing values :( \n {}'.format(np.array(atlas.labels[1:])[row_has_nan]))
ROI selected
 ['Left Insular Cortex' 'Right Insular Cortex'
 'Left Inferior Frontal Gyrus, pars triangularis'
 'Right Inferior Frontal Gyrus, pars triangularis'
 'Left Inferior Frontal Gyrus, pars opercularis'
 'Right Inferior Frontal Gyrus, pars opercularis'
 'Left Superior Temporal Gyrus, anterior division'
 'Right Superior Temporal Gyrus, anterior division'
 'Left Superior Temporal Gyrus, posterior division'
 'Right Superior Temporal Gyrus, posterior division'
 'Left Middle Temporal Gyrus, anterior division'
 'Right Middle Temporal Gyrus, anterior division'
 'Left Middle Temporal Gyrus, posterior division'
 'Right Middle Temporal Gyrus, posterior division'
 'Left Middle Temporal Gyrus, temporooccipital part'
 'Right Middle Temporal Gyrus, temporooccipital part'
 'Left Inferior Temporal Gyrus, temporooccipital part'
 'Right Inferior Temporal Gyrus, temporooccipital part'
 'Left Supramarginal Gyrus, anterior division'
 'Right Supramarginal Gyrus, anterior division' 'Left Angular Gyrus'
 'Left Lateral Occipital Cortex, inferior division'
 'Right Lateral Occipital Cortex, inferior division'
 'Left Intracalcarine Cortex' 'Right Intracalcarine Cortex'
 'Left Frontal Medial Cortex' 'Right Frontal Medial Cortex'
 'Left Subcallosal Cortex' 'Right Subcallosal Cortex'
 'Left Paracingulate Gyrus' 'Right Paracingulate Gyrus'
 'Left Cingulate Gyrus, anterior division'
 'Right Cingulate Gyrus, anterior division'
 'Left Cingulate Gyrus, posterior division'
 'Right Cingulate Gyrus, posterior division' 'Left Cuneal Cortex'
 'Right Cuneal Cortex' 'Left Frontal Orbital Cortex'
 'Right Frontal Orbital Cortex'
 'Left Parahippocampal Gyrus, anterior division'
 'Right Parahippocampal Gyrus, anterior division'
 'Left Parahippocampal Gyrus, posterior division'
 'Right Parahippocampal Gyrus, posterior division' 'Left Lingual Gyrus'
 'Right Lingual Gyrus'
 'Right Temporal Fusiform Cortex, posterior division'
 'Left Temporal Occipital Fusiform Cortex'
 'Right Temporal Occipital Fusiform Cortex'
 'Left Occipital Fusiform Gyrus' 'Right Occipital Fusiform Gyrus'
 'Left Frontal Operculum Cortex' 'Right Frontal Operculum Cortex'
 'Left Central Opercular Cortex' 'Right Central Opercular Cortex'
 'Left Parietal Operculum Cortex' 'Right Parietal Operculum Cortex'
 'Left Planum Polare' 'Right Planum Polare'
 "Left Heschl's Gyrus (includes H1 and H2)"
 "Right Heschl's Gyrus (includes H1 and H2)" 'Left Planum Temporale'
 'Right Planum Temporale' 'Left Supracalcarine Cortex'
 'Right Supracalcarine Cortex']
ROI removed due to missing values :( 
 ['Left Frontal Pole' 'Right Frontal Pole' 'Left Superior Frontal Gyrus'
 'Right Superior Frontal Gyrus' 'Left Middle Frontal Gyrus'
 'Right Middle Frontal Gyrus' 'Left Precentral Gyrus'
 'Right Precentral Gyrus' 'Left Temporal Pole' 'Right Temporal Pole'
 'Left Inferior Temporal Gyrus, anterior division'
 'Right Inferior Temporal Gyrus, anterior division'
 'Left Inferior Temporal Gyrus, posterior division'
 'Right Inferior Temporal Gyrus, posterior division'
 'Left Postcentral Gyrus' 'Right Postcentral Gyrus'
 'Left Superior Parietal Lobule' 'Right Superior Parietal Lobule'
 'Left Supramarginal Gyrus, posterior division'
 'Right Supramarginal Gyrus, posterior division' 'Right Angular Gyrus'
 'Left Lateral Occipital Cortex, superior division'
 'Right Lateral Occipital Cortex, superior division'
 'Left Juxtapositional Lobule Cortex (formerly Supplementary Motor Cortex)'
 'Right Juxtapositional Lobule Cortex (formerly Supplementary Motor Cortex)'
 'Left Precuneous Cortex' 'Right Precuneous Cortex'
 'Left Temporal Fusiform Cortex, anterior division'
 'Right Temporal Fusiform Cortex, anterior division'
 'Left Temporal Fusiform Cortex, posterior division' 'Left Occipital Pole'
 'Right Occipital Pole']

3.2 Compute FC and ISFC

Here we compute FC and ISFC on the parcellated data.

In [17]:
# Compute FC
fc_maps = {
    task_name:np.zeros((n_roi_select,n_roi_select)) 
    for task_name in all_task_names
}
for task_name in all_task_names: 
    for subj_id in range(n_subjs[task_name]): 
        fc_maps[task_name] += np.corrcoef(
            bold_ho_filtered[task_name][:,:,subj_id].T
        ) / n_subjs[task_name]
        np.fill_diagonal(fc_maps[task_name], np.nan)     
In [18]:
# Compute ISFC
isfc_maps_ho = {}
for task_name in all_task_names:
    isfc_maps_ho[task_name] = isfc(data=bold_ho_filtered[task_name],
                                   summary_statistic='median',
                                   vectorize_isfcs=False)   

Exercise 6: Visualize the FC/ISFC matrices, averaged across subjects

  • Use imshow to visualize the 4 correlation matrices: (FC vs. ISFC) x (word-level scrambled vs. intact story)
  • Mark the rows (or columns) with the ROI labels
In [19]:
# Insert code here

Exercise 7: Visualize FC/ISFC connectivity strength on glass brains

  • Use plot_connectome to visualize the 4 correlation matrices: (FC vs. ISFC) x (word-level scrabled vs. intact story)
  • Use common edge_threshold for all plots.

Hint: The plot_connectome function takes as an input a correlation matrix (such as the FC one plotted above) and also a set of coordinates that define the XYZ coordinates that corresponds to each column/row of the correlation matrix. In order to get the coordinates of these ROIs, we recommend you use: plotting.find_parcellation_cut_coords(atlas.maps)[roi_select]

In [20]:
# Insert code here

Exercise 8: Do FC maps look different across conditions? How about ISFC? And why? Hint: consult this paper.

A:

4. Spatial pattern correlation across subjects

4.1 Spatial inter-subject correlation


We can apply the idea of inter-subject analysis to RSA. So far, ISC is being computed between aligned pairs of voxels across time points and it is commonly referred to as temporal ISC. However, we could instead correlate between aligned pairs of time points across voxels. That is, how does the pattern of activity across voxels for one time point correlate with the average pattern of the other participants at that time point. By doing this for each time point, we can generate a time course of these correlations to observe the general ebb and flow of coupling in brain activity across participants. This can be done simply by transposing the voxel and time dimensions (for a 3-D matrix, this is accomplished with a 90 degree rotation). This is a simple transposition. If we have data in the format: (TRs, voxels, subjects), we can use data.transpose(1,0,2), where the indices refer to the dimensions of the array and we will have an array in the format (voxels, TRs, subjects).

Compare the two task conditions with spatial ISC

One way to compare the intact and word_scramble conditions is to plot the correlation values, by TR, for each condition. Again, we can use the same ISC functions from above after transposing the data matrices.

In [21]:
# Get a list of ROIs. 
roi_mask_path = os.path.join(pieman2_dir,'masks','rois')
all_roi_fpaths = glob.glob(os.path.join(roi_mask_path, '*.nii.gz'))

# Collect all ROIs 
all_roi_names = []
all_roi_nii = {}
all_roi_masker = {}
for roi_fpath in all_roi_fpaths:
    
    # Compute ROI name
    roi_fname = os.path.basename(roi_fpath)
    roi_name = roi_fname.split('.')[0]
    all_roi_names.append(roi_name)
    
    # Load roi nii file 
    roi_nii = nib.load(roi_fpath)
    all_roi_nii[roi_name] = roi_nii
    
    # Make roi maskers
    all_roi_masker[roi_name] = NiftiMasker(mask_img=roi_nii)

print('Path to all roi masks: {}'.format(roi_mask_path))    
print('Here are all ROIs:\n{}'.format(all_roi_names))
Path to all roi masks: /jukebox/norman/mkumar/tutorial_data/brainiak_datasets/Pieman2/masks/rois
Here are all ROIs:
['aMTG_R', 'aIFG_R', 'MCC', 'caMTG_L', 'hV4_R', 'SOG_L', 'dMFG_R', 'lOFG_R', 'aOFC_R', 'Insular_L', 'SMG_R', 'vmPFC', 'dmPFG', 'dPoCS_L', 'V3a_L', 'vPCUN', 'aINS_L', 'dMFG_L', 'pMTG_R', 'lOFG_L', 'V1', 'ldV2_L', 'aPFC', 'IPL_L', 'STC_R', 'DLPFC_L', 'V3a_R', 'DLPFC_R', 'aANG_L', 'vPFC_R', 'LOC_L', 'pANG_L', 'Insular_R', 'aSTG_L', 'smPFC', 'dPCUN', 'aIFG_L', 'SMA_L', 'STC_L', 'IPL_R', 'HG_L', 'pMTS_L', 'aANG_R', 'aCUN', 'V3', 'pMTS_R', 'LOC_R', 'mldV2_R', 'pANG_R', 'vV2', 'SFG', 'cMTG_R', 'SFS_L', 'vMFG_L', 'dPoCS_R', 'pIFG_L', 'hV4_L', 'PMC_L', 'HG_R', 'dPCC', 'aOFC_L', 'SMG_L', 'dPreCG_L']
In [22]:
# Make a function to load data for one ROI
def load_roi_data(roi_name): 
    # Pick a roi masker
    roi_masker = all_roi_masker[roi_name]    
    
    # Preallocate 
    bold_roi = {task_name:[] for i, task_name in enumerate(all_task_names)}
    
    # Gather data 
    for task_name in all_task_names:
        for subj_id in range(n_subjs[task_name]):
            
            # Get the data for task t, subject s 
            nii_t_s = nib.load(fnames[task_name][subj_id])
            bold_roi[task_name].append(roi_masker.fit_transform(nii_t_s))
            
        # Reformat the data to std form 
        bold_roi[task_name] = np.transpose(np.array(bold_roi[task_name]), [1,2,0])
    return bold_roi

Compute spatial ISC on some ROIs.

In [23]:
roi_selected = ['dPCC', 'vPCUN', 'V1']
roi_selected_names = ['dorsal posterior cingulate cortex', 'ventral precuneus', 'primary visual cortex']

# compute sISC for all ROIs 
iscs_roi_selected = []
for j, roi_name in enumerate(roi_selected):
    print(j, roi_name)
    
    # Load data 
    bold_roi = load_roi_data(roi_name)
    
    # Compute isc 
    iscs_roi = {}
    for task_name in all_task_names: 
        iscs_roi[task_name] = isc(np.transpose(bold_roi[task_name], [1,0,2]))
        
    iscs_roi_selected.append(iscs_roi)
0 dPCC
1 vPCUN
2 V1
In [24]:
# Plot the spatial ISC over time
col_pal = sns.color_palette(palette='colorblind', n_colors=len(all_task_names))
ci = 95

f, axes = plt.subplots(len(roi_selected), 1, figsize=(14, 5 * len(roi_selected)), sharex=True)

# For each ROI
for j, roi_name in enumerate(roi_selected):
    # For each task 
    for i, task_name in enumerate(all_task_names): 
        sns.tsplot(
            iscs_roi_selected[j][task_name], 
            color=col_pal[i], ci=ci, 
            ax=axes[j]
        )
    f.legend(all_task_des)
    sns.despine()

# Label the plot 
for j, roi_name in enumerate(roi_selected):
    axes[j].axhline(0, color='black', linestyle='--', alpha=.3)
    axes[j].set_ylabel('Linear correlation')
    axes[j].set_title('Spatial inter-subject correlation, {}'. format(roi_selected_names[j]))
    
axes[-1].set_xlabel('TRs')
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/usr/people/mk35/.conda/envs/mybrainiak/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
Out[24]:
Text(0.5, 0, 'TRs')

Exercise 9: Interpret the spatial ISC results you observed above.

A:

Novel contribution: be creative and make one new discovery by adding an analysis, visualization, or optimization.

Here are some ideas:

  • The ISC package in BrainIAK supports other statistical tests. Study one of them, describe the logic behind it, and re-run the analysis with that test.
  • Conduct a sliding window spatial ISC analysis.
  • Perform some clustering on the ISFC matrices. Use the clustering (instead of the anatomical parcel) to re-analyze the data, and compare ISFC vs. FC across tasks. See here for a good resource on clustering algorithms.

Contributions

E. Simony and U. Hasson for providing data
C. Baldassano and C. Chen provided initial code
M. Kumar, C. Ellis and N. Turk-Browne produced the initial notebook 4/4/18
S. Nastase enhanced the ISC brainiak module; added the section on statistical testing
Q. Lu added solutions; switched to S. Nastase's ISC module; replicated Lerner et al 2011 & Simony et al. 2016; added spatial ISC.
M. Kumar edits to section introductions and explanation on permutation test.
K.A. Norman provided suggestions on the overall content and made edits to this notebook.
C. Ellis incorporated edits from cmhn-s19