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

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")
<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.)')

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 = (, 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)")
<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')
<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)
<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
# 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")
<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!
 <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))
<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): 
# 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")
<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')
<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')
<nilearn.plotting.displays.OrthoSlicer at 0x7fec81e63ed0>