Note

Go to the end to download the full example code. or to run this example in your browser via Binder

Regions extraction using dictionary learning and functional connectomes¶

This example shows how to use RegionExtractor to extract spatially constrained brain regions from whole brain maps decomposed using Dictionary learning and use them to build a functional connectome.

We used movie-watching functional scans of 20 subjects from fetch_development_fmri and DictLearning for set of brain atlas maps.

This example can also be inspired to apply the same steps to even regions extraction using ICA maps. In that case, the idea would be to replace Dictionary learning to canonical ICA decomposition using CanICA

Please see the related documentation of RegionExtractor for more details.

Fetch brain development functional datasets¶

We use nilearn’s dataset downloading utilities:

from nilearn.datasets import fetch_development_fmri

rest_dataset = fetch_development_fmri(n_subjects=20)
func_filenames = rest_dataset.func
confounds = rest_dataset.confounds

[fetch_development_fmri] Dataset directory found: /home/runner/work/nilearn/nilearn/nilearn_data/development_fmri

Extract functional networks with Dictionary learning ¶

Import DictLearning from the decomposition module, instantiate the object, and fit the model to the functional datasets

from nilearn.decomposition import DictLearning

# Initialize DictLearning object
dict_learn = DictLearning(
    n_components=8,
    smoothing_fwhm=6.0,
    memory="nilearn_cache",
    memory_level=1,
    random_state=0,
    verbose=1,
)
# Fit to the data
dict_learn.fit(func_filenames)

[DictLearning.fit] Loading data
[DictLearning.fit] Learning initial components
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    1.3s finished
[DictLearning.fit] Computing initial loadings
________________________________________________________________________________
[Memory] Calling nilearn.decomposition.dict_learning._compute_loadings...
_compute_loadings(array([[-0.006357, ..., -0.000176],
       ...,
       [ 0.005972, ..., -0.001489]]),
array([[-3.895898, ..., -0.034474],
       ...,
       [ 1.466817, ..., -0.956513]]))
_________________________________________________compute_loadings - 0.0s, 0.0min
[DictLearning.fit]  Learning dictionary
________________________________________________________________________________
[Memory] Calling sklearn.decomposition._dict_learning.dict_learning_online...
dict_learning_online(array([[-3.895898, ...,  1.466817],
       ...,
       [-0.034474, ..., -0.956513]]),
8, alpha=10, batch_size=20, method='cd', dict_init=array([[ 0.288595, ..., -0.536917],
       ...,
       [-0.154675, ..., -0.35914 ]]), verbose=0, random_state=0, return_code=True, shuffle=True, n_jobs=1, max_iter=1414)
_____________________________________________dict_learning_online - 0.9s, 0.0min

DictLearning(memory='nilearn_cache', memory_level=1, n_components=8,
             random_state=0, smoothing_fwhm=6.0, verbose=1)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Visualization of functional networks¶

Get the extracted functional networks via the attribute components_img_ and show them using plotting utilities.

from nilearn.plotting import plot_prob_atlas, show

components_img = dict_learn.components_img_

plot_prob_atlas(
    components_img,
    view_type="filled_contours",
    title="Dictionary Learning maps",
    draw_cross=False,
)

show()

Extract regions from networks¶

Import RegionExtractor from the regions module. threshold=0.5 indicates that we keep nominal of amount nonzero voxels across all maps, less the threshold means that Import RegionExtractor from the regions module. threshold=0.5 indicates that we keep a nominal of amount non-zero voxels across all maps.

from nilearn.regions import RegionExtractor

extractor = RegionExtractor(
    components_img,
    threshold=0.5,
    standardize="zscore_sample",
    thresholding_strategy="ratio_n_voxels",
    extractor="local_regions",
    standardize_confounds=True,
    min_region_size=1350,
    verbose=1,
)

# Just call fit() to proceed with regions extraction.
extractor.fit()

[RegionExtractor.fit] Loading regions from <nibabel.nifti1.Nifti1Image object at 0x7f17bcdebd50>
[RegionExtractor.fit] Finished fit

RegionExtractor(maps_img=<nibabel.nifti1.Nifti1Image object at 0x7f1804d5c150>,
                standardize='zscore_sample', threshold=0.5, verbose=1)

Visualization of region extraction results¶

Extracted regions are stored in regions_img_ and each region index is stored in index_.

regions_index = extractor.index_

regions_extracted_img = extractor.regions_img_
n_regions_extracted = regions_extracted_img.shape[-1]
n_components = components_img.shape[-1]
title = (
    f"{n_regions_extracted} regions are extracted "
    f"from {n_components} components.\n"
    "Each separate color of region indicates extracted region."
)

plot_prob_atlas(
    regions_extracted_img,
    view_type="filled_contours",
    title=title,
    draw_cross=False,
)

show()

/home/runner/work/nilearn/nilearn/.tox/doc/lib/python3.11/site-packages/numpy/ma/core.py:2820: UserWarning:

Warning: converting a masked element to nan.

Compute correlation coefficients¶

First we need to do subjects timeseries signals extraction and then estimating correlation matrices on those signals. To extract timeseries signals, we call transform onto each subject functional data stored in func_filenames. First we need to extract timeseries signals for each subject and then estimate the correlation matrices on those signals. To extract timeseries, we call transform onto each subject functional data stored in func_filenames. To estimate the correlation matrices we can rely on nilearn.connectome.ConnectivityMeasure.

from nilearn.connectome import ConnectivityMeasure

connectome_measure = ConnectivityMeasure(kind="correlation", verbose=1)

correlations = []
for filename, confound in zip(func_filenames, confounds, strict=False):
    timeseries_each_subject = extractor.transform(filename, confounds=confound)
    correlation = connectome_measure.fit_transform([timeseries_each_subject])
    correlations.append(correlation)

[RegionExtractor.wrapped] Loading data from sub-pixar126_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar124_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar123_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar125_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar016_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar015_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar014_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar013_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar012_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar011_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar001_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar008_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar007_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar006_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar005_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar004_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar003_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar002_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar009_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit
[RegionExtractor.wrapped] Loading data from sub-pixar010_task-...
[RegionExtractor.wrapped] Smoothing images
[RegionExtractor.wrapped] Extracting region signals
[RegionExtractor.wrapped] Cleaning extracted signals
[ConnectivityMeasure.wrapped] Finished fit

Plot resulting connectomes¶

First we plot the mean of correlation matrices with plot_matrix, and we use plot_connectome to plot the connectome relations.

import numpy as np

from nilearn.plotting import (
    find_probabilistic_atlas_cut_coords,
    find_xyz_cut_coords,
    plot_connectome,
    plot_matrix,
)

mean_correlations = np.mean(correlations, axis=0).reshape(
    n_regions_extracted, n_regions_extracted
)

title = f"Correlation between {int(n_regions_extracted)} regions"
plot_matrix(mean_correlations, vmax=1, vmin=-1, title=title)

# Find the center of the regions and plot the connectome.
regions_img = regions_extracted_img
coords_connectome = find_probabilistic_atlas_cut_coords(regions_img)
plot_connectome(
    mean_correlations, coords_connectome, edge_threshold="90%", title=title
)

show()

Plot regions extracted for only one specific network¶

First, we plot a network of index=4 without region extraction.

from nilearn import image
from nilearn.plotting import plot_stat_map

index = 4

img = image.index_img(components_img, index)
coords = find_xyz_cut_coords(img)
plot_stat_map(
    img,
    cut_coords=coords,
    title="Showing one specific network",
)

show()

Now, we plot the same network after region extraction to show that connected regions are nicely separated. Each brain extracted region is identified as separate color.

For this, we take the indices of the all regions extracted related to original network given as 4.

from nilearn.plotting import cm, plot_anat

regions_indices_of_map3 = np.where(np.array(regions_index) == index)

display = plot_anat(
    cut_coords=coords, title="Regions from this network", colorbar=False
)

# Add as an overlay all the regions of index 4
colors = "rgbcmyk"
for each_index_of_map3, color in zip(
    regions_indices_of_map3[0], colors, strict=False
):
    display.add_overlay(
        image.index_img(regions_extracted_img, each_index_of_map3),
        cmap=cm.alpha_cmap(color),
    )

show()

Total running time of the script: (0 minutes 59.140 seconds)

Estimated memory usage: 1285 MB

Gallery generated by Sphinx-Gallery

Regions extraction using dictionary learning and functional connectomes¶

Fetch brain development functional datasets¶

Extract functional networks with Dictionary learning¶

Visualization of functional networks¶

Extract regions from networks¶

Visualization of region extraction results¶

Compute correlation coefficients¶

Plot resulting connectomes¶

Plot regions extracted for only one specific network¶

Extract functional networks with Dictionary learning ¶