Introduction to Functional Brain Parcellation

Pierre Bellec

Why parcellations? understand brain organization

Local and distributed connectivity lead to the emergence of both local and distributed neuronal assemblies, also called functional networks. Figures from Varela et al., 2001.

Why parcellations? reduce dimensionality

Parcels can be used as nodes to approximate brain networks as graphs. Here, the average connectome of the ADHD-200 sample is represented with several different parcellations, of varying resolutions. From Bellec et al., Neuroimage (2017).

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

Vasculature meshes with neuronal population at a very fine spatial scale, (~10 microns), with micro-capilaries regulating blood oxygenation in a highly local and precise way. The Blood Oxygenation-Level Dependent (BOLD) signal in fMRI captures this neurovasculare coupling. Figure adapted from Harrison, 2002 and Dr Bruce Pike. Check this lecture for more info.

functional Magnetic Resonance Imaging

fMRI is a 4D imaging modality. Each brain parcel, down to a single volume element (voxel), is associated with a time series.

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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec8a0d2310>

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))")
Time series with shape (176, 444) (# 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

Slow spontaneous fluctuations and seed-based connectivity map from the posterior cingulate cortex identifies the default-mode network. See the following nilearn tutorial for more details.

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))")
Time series with shape (176, 69681) (# 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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec850bc4d0>

Functional Connectome

A functional connectome essentially is a collection of seed-based connectivity maps, approximated (or compressed) on parcels, and using all possible seeds.

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]:
<matplotlib.image.AxesImage at 0x7fec84a98c50>

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]:
<nilearn.plotting.displays.XSlicer at 0x7fec83f9ba50>

Functional Networks

With the rigth order on parcels, we can see functional networks as diagonal squares with high connectivity on the diagonal. Finding these squares (and the order) can be solved with a cluster analysis.

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
[372, 385, 4, 53, 197, 287, 313, 346, 62, 192]

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?
[ 0  1  2  3  4  5  5  6  7  8  8  5  9  0 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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec81eaff50>

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]:
(<AxesSubplot:title={'center':'part'}>,
 <matplotlib.image.AxesImage at 0x7fec82377e50>,
 Text(0.5, 1.0, 'part'))

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]:
<nilearn.plotting.displays.XSlicer at 0x7fec81ee8bd0>

From individual to group parcels

Adapted from Liu et al., Neuroimage (2014).

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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec8212c0d0>

Multiresolution parcellations

Golland parcellation (2 clusters)

Networks can also be merged into reliable “exogeneous” vs “endogeneous” systems. From Golland et al. (2008).

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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec822e04d0>

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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec81e63ed0>

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)

Homogeneity of group brain parcels as a function of the average parcel size. Note that homogeneity can closely be predicted from size alone. From Urchs et al. MNI open research (2019).

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]:
<nilearn.plotting.displays.XSlicer at 0x7fec81ab1150>

Group reproducibility

With large, homogeneous groups, group-level parcellations reproduce closely. Figure adapted from Yeo, Krienen et al. (2011).

Individual parcels

For ten densely sampled individuals (10 runs of 30 mns resting-state over ten days), Gordon et al. (2017) identified details in individual parcellations that cannot be observed at the level of group parcellations (indicated by arrows, group parcellation at the top).

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)
adjusted Rand reproducibility: 0.1058196841391461

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

Figure from the excellent introduction to multivariate analysis of fMRI by (Khosla et al., ArXiv 2018)

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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec81b5bbd0>
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]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fec819eec10>

Which one is best?

Sparse dictionary learning has a (short) edge in classification benchmarks.

Figure from (Dadi et al., ArXiv 2020)

Which one is best?

Sparse dictionary learning has a (clear) edge in compression quality.

Figure from (Dadi et al., ArXiv 2020)

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.

Keep calm and parcellate!

Art by Emilio Garcia

Powered by nilearn