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]:
Multiresolution parcellations¶
Some parcellation atlases explicitely offer a range of resolution, such as the MIST atlas below (Urchs et al., 2019) (see dashboard) and the Schaefer atlas (Schaefer et al., 2018).
In [19]:
plt.figure(figsize=(18, 6))
for i, scale in enumerate(basc):
if not scale == 'description':
plotting.plot_stat_map(basc[scale], display_mode="x", title=scale, cut_coords=(0, ),
colorbar=False, cmap="nipy_spectral", axes=plt.subplot(2, 5, i + 1))
Cluster homogeneity (compression quality)
In [20]:
# Let's check qualitatively how compression quality varies with scale.
plt.figure(figsize=(18, 4))
for i, scale in enumerate(("scale020", "scale064", "scale444")):
masker_scale = input_data.NiftiLabelsMasker(basc[scale], resampling_target="data", detrend=True,
standardize=True, smoothing_fwhm=5, memory='nilearn_cache', memory_level=1).fit()
tseries_scale = masker_scale.transform(adhd.func[0])
seed_to_voxel_correlations = (np.dot(tseries_scale.T, tseries[:, 316]) / tseries.shape[0])
conn_map_scale = masker_scale.inverse_transform(seed_to_voxel_correlations.reshape([1, tseries_scale.shape[1]]))
plotting.plot_stat_map(conn_map_scale, display_mode='x', cut_coords=(2, ), threshold=0.5, vmax=1,
title=f"{scale} fcMRI map (PCC)", axes=plt.subplot(1, 4, i + 1), colorbar=False)
plotting.plot_stat_map(conn_map, display_mode='x', cut_coords=(2, ), threshold=0.5, vmax=1,
title="voxel fcMRI map (PCC)", axes=plt.subplot(1, 4, 4), colorbar=False)
Out[20]:
Group reproducibility
With large, homogeneous groups, group-level parcellations reproduce closely. Figure adapted from Yeo, Krienen et al. (2011).
Individual parcels
Individual reproducibility
Even for long time series, the reproducibility of individual parcellations seems to be reaching a plateau (Gordon et al., 2017).
In [21]:
# Let's compare two parcellations generated from different subjects
part_ind = []
for ind in range(2):
tseries_ind = masker.transform(adhd.func[ind])
conn_ind = np.squeeze(ConnectivityMeasure(kind='correlation').fit_transform([tseries_ind]))
hier_ind = linkage(conn_ind, method='average', metric='euclidean')
part_ind.append(np.squeeze(cut_tree(hier_ind, n_clusters=17)))
# If we were to count the proportion of elements in the two adjacency matrices that are identical
# (excluding the diagonal which is always 1), we would get a measure of agreement called the Rand index.
# The adjusted Rand score corrects for chance-level overlap, such that the adjusted Rand will be close to
# zero for cluster overlap near chance level, and adjusted Rand is 1 for identical cluster solutions.
# Note that scikit-learn actually does not need the adjacency matrix representation, but works directly
# from vectors of cluster labels.
from sklearn import metrics
repro = metrics.adjusted_rand_score(part_ind[0], part_ind[1])
print("adjusted Rand reproducibility:", repro)
Ouch, this is, like, really low. Inter-subject variability? Not enough data? Bad algorithm? all of the above? Check (Nikolaidis et al., Neuroimage 2020) for a systematic evaluation.
Linear decomposition models
In [22]:
# Let's use dictionary learning from nilearn
from nilearn.decomposition import DictLearning
# note that the API is quite standard, follows sklearn
# Pretty much all the methods covered so far can be run that easily!
# Note this interface includes the options we previously used for the masker
dict_learn = DictLearning(n_components=17, detrend=True, standardize=False, smoothing_fwhm=5,
memory="nilearn_cache", memory_level=2, random_state=0)
dict_learn.fit(adhd.func[0])
# Let's have a look with the prob atlas viewer
plotting.plot_prob_atlas(dict_learn.components_img_, view_type='filled_contours', title='Dictionary Learning maps')
Out[22]:
In [23]:
# Now let's compare the DMN cluster with the DMN component
from nilearn import image
plotting.plot_stat_map(image.index_img(dict_learn.components_img_, 13), threshold=0.01, title="DMN component",
cut_coords=(2, -51, 27))
plotting.plot_roi(part_img, title="DMN cluster", cut_coords=(2, -51, 27))
# It does not always work that well!! There are qualitatively different networks that component methods
# and cluster analysis identify
Out[23]:
Which one is best?
Sparse dictionary learning has a (short) edge in classification benchmarks.
Which one is best?
Sparse dictionary learning has a (clear) edge in compression quality.
Summary¶
- Functional parcelations are modules of brain regions with higher homogeneity intra-module than inter-modules.
- Parcelations exist over a range of scales, from distributed network down to specialized cortical areas.
- Homogeneity increases with scale.
- Stability is high at group-level, uneven at individual-level.
- A number of benchmarks are now established to compare parcellation algorithms.
- There are substantial differences between individual and group parcellations, and these two levels can be estimated jointly.
- Gradients (or any linear mixture) provide a complementary, richer view on brain organization than parcels, but are also more complex.
Powered by nilearn
