Introduction to Functional Brain Parcellation
Why parcellations? understand brain organization
Why parcellations? reduce dimensionality
Code and data¶
Nilearn code is presented to generate figures similar to most slides in this presentation. To run the code locally, clone the github repository and install the dependencies listed in requirements.txt and run the following code to cache all the necessary data:
In [1]:
import warnings
warnings.filterwarnings("ignore")
from nilearn import datasets # Fetch data using nilearn
basc = datasets.fetch_atlas_basc_multiscale_2015() # the BASC multiscale atlas
adhd = datasets.fetch_adhd(n_subjects=10) # ADHD200 preprocessed data (Athena pipeline)
atlas_yeo = datasets.fetch_atlas_yeo_2011() # the Yeo-Krienen atlas
The BOLD signal
functional Magnetic Resonance Imaging
Brain parcels¶
In [2]:
# Let's load some small brain parcels!
from nilearn import plotting, input_data
plotting.plot_roi(basc["scale444"], title="BASC scale 444", colorbar=True, cmap="nipy_spectral")
Out[2]:
Parcel time series extraction¶
In [3]:
# We use a nilearn masker to load time series from the parcellation above
masker = input_data.NiftiLabelsMasker(basc['scale444'], resampling_target="data", detrend=True,
standardize=True, smoothing_fwhm=5, memory='nilearn_cache', memory_level=1).fit()
tseries = masker.transform(adhd.func[0])
print(f"Time series with shape {tseries.shape} (# time points, # parcels))")
In [4]:
# plot time series for one region
import matplotlib.pyplot as plt
plt.plot(tseries[:,316]), plt.title('Seed time series (PCC)'), plt.xlabel('Time'), plt.ylabel('BOLD (a.u.)')
plt.tight_layout()
Seed-based connectivity map
Loading voxel data¶
In [5]:
# Load voxel-level data - a brain mask is automatically computed
masker_voxel = input_data.NiftiMasker(detrend=True, standardize=True, smoothing_fwhm=5).fit(adhd.func[0])
tseries_voxel = masker_voxel.transform(adhd.func[0])
print(f"Time series with shape {tseries_voxel.shape} (# time points, # voxels))")
Computing the voxel-based connectivity map¶
In [6]:
# Compute correlation, and look at voxel-based connectivity map
import numpy as np
seed_to_voxel_correlations = (np.dot(tseries_voxel.T, tseries[:, 316]) / tseries.shape[0])# Show the connectivity map
conn_map = masker_voxel.inverse_transform(seed_to_voxel_correlations.T)
plotting.plot_stat_map(conn_map, cut_coords=(2, -51, 27), threshold=0.5, vmax=1, title="fcMRI map (PCC)")
Out[6]:
Functional Connectome
Functional Connectome
In [7]:
# We generate a connectome
from nilearn.connectome import ConnectivityMeasure
conn = np.squeeze(ConnectivityMeasure(kind='correlation').fit_transform([tseries]))
# Plot the connectome, with no particular region order
plotting.plot_matrix(conn, labels = range(444), title='Functional connectome', reorder='average')
Out[7]:
Compressed connectivity map¶
(Craddock et al., 2012) proposed to evaluate parcellation through compression.
In [8]:
# Let's project back one of the column of the connectome in voxel space, and show the original voxel-based map for comparison
conn_map_compressed = masker.inverse_transform(conn[316, :].reshape([1, 444]))
plotting.plot_stat_map(conn_map_compressed, display_mode='x', cut_coords=(2, ), axes=plt.subplot(1, 2, 1), threshold=0.5, vmax=1, title="compressed fcMRI map (PCC)", colorbar=False)
plotting.plot_stat_map(conn_map, display_mode='x', cut_coords=(2, ), axes=plt.subplot(1, 2, 2), threshold=0.5, vmax=1, title="fcMRI map (PCC)", colorbar=False)
Out[8]:
Functional Networks
Hierarchical cluster analysis: growing a tree¶
In [9]:
# we use scipy's hierarchical clustering implementation
from scipy.cluster.hierarchy import dendrogram, linkage, cut_tree
# That's the hierarchical clustering step
hier = linkage(conn, method='average', metric='euclidean') # scipy's hierarchical clustering
# HAC proceeds by iteratively merging brain regions, which can be visualized with a tree
res = dendrogram(hier, get_leaves=True) # Generate a dendrogram from the hierarchy
Hierarchical cluster analysis: ordering leaves¶
In [10]:
# the order of merging above give us a good order to visualize the matrix
order = res.get('leaves') # Extract the order on parcels from the dendrogram
print(order[0:10])
# here parcel 372 is the left-most leaf in the tree, then 385, etc
Hierarchical cluster analysis: cutting a tree¶
In [11]:
# We can cut the tree at whatever number of clusters we choose (here 17, because 17 is good)
part = np.squeeze(cut_tree(hier, n_clusters=17)) # Cut the hierarchy
# Each entry of the vector part is a parcel, and codes for the number of the network of this parcel
print(part[0:15]) # e.g. parcel #7 is in cluster 5. What is the cluster of parcel number 10?
In [12]:
# let's visualize these clusters back in voxel space
part_img = masker.inverse_transform(part.reshape([1, 444]) + 1) # note the sneaky shift to 1-indexing
plotting.plot_roi(part_img, title="functional networks", colorbar=True, cmap="Paired")
Out[12]:
Hierarchical cluster analysis: order and adjacency¶
In [13]:
# Another way to represent the clusters is an adjacency matrix A
# A(i,j) = 1 if i and j are in the same cluster, and 0 otherwise
def part2adj(part):
part = np.reshape(part, [part.shape[0], 1])
adj = np.repeat(part, part.shape[0], 1) == np.repeat(part.transpose(), part.shape[0], 0)
return adj.astype(int)
In [14]:
# Let's look at the adjacency matrix
# alongside the connectivity matrix, after reordering
plt.subplot(1,2,1), plt.imshow(part2adj(part[order])), plt.title('part')
plt.subplot(1,2,2), plt.imshow(conn[order, :][:, order], cmap='bwr'), plt.title('part')
# in the adjacency matrix, the clusters are literally squares on the diagonal!
Out[14]:
Hierarchical cluster analysis: compression¶
In [15]:
# Let's finally visualize the compressed seed-based PCC map
conn_map_compressed = masker.inverse_transform(conn[316, :].reshape([1, 444]))
plotting.plot_stat_map(conn_map_compressed, display_mode='x', cut_coords=(2, ), threshold=0.5, vmax=1,
axes=plt.subplot(1, 2, 1), title="PCC fcMRI (BASC444)", colorbar=False)
# along with the DMN cluster
part_img = masker.inverse_transform(part.reshape([1, 444]) + 1 == 10) # note the sneaky shift to 1-indexing
plotting.plot_roi(part_img, title="DMN cluster", display_mode='x', cut_coords=(2, ), axes=plt.subplot(1, 2, 2))
Out[15]:
In [16]:
# Let's do a crude group analysis: load time series for a few subjects (N=10)
tseries_all = []
for indx, img in enumerate(adhd.func):
tseries_all.append(masker.transform(img))
# clustering on the average connectome
conn_all = ConnectivityMeasure(kind='correlation').fit(tseries_all)
hier = linkage(conn_all.mean_, method='average', metric='euclidean')
# Cut the hierarchy and look at the group-level networks
part_group = np.squeeze(cut_tree(hier, n_clusters=17))
part_group_img = masker.inverse_transform(part_group.reshape([1, 444]) + 1) # sneaky shift to 1-indexing
plotting.plot_roi(part_group_img, title="group networks", colorbar=True, cmap="Paired")
Out[16]:
Multiresolution parcellations
Golland parcellation (2 clusters)
Yeo-Krienen parcellation (7 clusters)¶
Group cluster analysis of fMRI connectivity (N=1000) using von Mises-Fisher distribution (Lashkari et al., 2010). See Yeo, Krienen et al. (2011) and this tutorial. Also check as an interactive view.
In [17]:
# Let's plot the Yeo-Krienen 7 clusters parcellation
from nilearn import plotting
plotting.plot_roi(atlas_yeo.thick_7, title='Yeo-Krienen atlas scale 7',
cut_coords=(8, -4, 9), colorbar=True, cmap='Paired')
Out[17]:
Yeo-Krienen parcellation (17 clusters)¶
Networks can be further divided into subnetworks, e.g. in the visual cortex. From Yeo, Krienen et al. (2011). Also check as an interactive view.
In [18]:
# Let's plot the Yeo-Krienen 7 clusters parcellation
from nilearn import plotting
plotting.plot_roi(atlas_yeo.thick_17, title='Yeo-Krienen atlas scale 17',
cut_coords=(8, -4, 9), colorbar=True, cmap='Paired')
Out[18]: