Connectivity
Written by Luke Chang So far, we have primarily been focusing on analyses related to task evoked brain activity. However, an entirely different way to study the brain is to characterize how it is intrinsically connected. There are many different ways to study functional connectivity. The primary division is studying how brain regions are structurally connected. In animal studies this might involve directly tracing bundles of neurons that are connected to other neurons. Diffusion imaging is a common way in which we can map how bundles of white matter are connected to each region, based on the direction in which water diffuses along white matter tracks. There are many different techniques such as fractional ansiotropy and probablistic tractography. We will not be discussing structural connectivity in this course. An alternative approach to studying connectivity is to examine how brain regions covary with each other in time. This is referred to as functional connectivity, but it is better to think about it as temporal covariation between regions as this does not necessarily imply that two regions are directly communication with each other.
For example, regions can directly influence each other, or they can indirectly influence each other via a mediating region, or they can be affected similarly by a shared influence. These types of figures are often called graphs. These types of graphical models can be directed or undirected. Directed graphs imply a causal relationship, where one region A directly influence another region B. Directed graphs or causal models are typically described as effective connectivity, while undirected graphs in which the relationship is presumed to be bidirectional are what we typically describe as functional connectivity.
In this tutorial, we will work through examples on:
- Seed-based functional connectivity
- Psychophysiological interactions
- Principal Components Analysis
- Graph Theory
Now, let's dive in a little bit deeper into the specific details of functional connectivity.
Functional Connectivity¶
Seed Voxel Correlations¶
One relatively simple way to calculate functional connectivity is to compute the temporal correlation between two regions of interest (ROIs). Typically, this is done by extracting the temporal response from a seed voxel or the average response within a seed region. Then this time course is regressed against all other voxels in the brain to produce a whole brain map of anywhere that shares a similar time course to the seed.
Let's try it ourselves with an example subject from the Pinel Localizer dataset. First, let's import the modules we need for this tutorial and set our paths.
# '%matplotlib inline' command supported automatically in marimo
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.fft import fft, fftfreq
from nltools.data import Brain_Data, Design_Matrix, Adjacency
from nltools.mask import expand_mask, roi_to_brain
from nltools.stats import zscore, fdr, one_sample_permutation
from nltools.file_reader import onsets_to_dm
from scipy.stats import binom, ttest_1samp
from sklearn.metrics import pairwise_distances
from copy import deepcopy
import networkx as nx
from nilearn.plotting import plot_stat_map, view_img_on_surf
import nibabel as nib
from dartbrains_tools.data import get_file, get_tr, load_events, load_confounds, get_subjects
from dartbrains_tools.notebook_utils import youtube
def get_csf_mask_path(subject):
"""Per-subject CSF probability mask from fmriprep outputs."""
from pathlib import Path as _Path
bold_path = _Path(get_file(subject, 'derivatives', 'bold'))
return str(bold_path.parent.parent / 'anat' /
f'sub-{subject}_space-MNI152NLin2009cAsym_label-CSF_probseg.nii.gz')
Now let's load an example participant's preprocessed functional data.
sub = 'S01'
_fwhm = 6
data = Brain_Data(get_file(sub, 'derivatives', 'bold'))
with mo.persistent_cache(name="connectivity_smoothed"):
smoothed = data.smooth(fwhm=_fwhm)
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:253: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = self.nifti_masker.fit_transform(data) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:2369: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. out.data = out.nifti_masker.fit_transform(smooth_img(self.to_nifti(), fwhm))
Next we need to pick an ROI. Pretty much any type of ROI will work.
In this example, we will be using a whole brain parcellation based on similar patterns of coactivation across over 10,000 published studies available in neurosynth (see this paper for more details). We will be using a parcellation of 50 different functionally similar ROIs.
with mo.persistent_cache(name="connectivity_k50_mask"):
mask_1 = Brain_Data('https://neurovault.org/media/images/8423/k50_2mm.nii.gz')
mask_1.iplot()
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:253: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = self.nifti_masker.fit_transform(data)
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order)
Each ROI in this parcellation has its own unique number. We can expand this so that each ROI becomes its own binary mask using nltools.mask.expand_mask.
Let's plot the first 5 masks.
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:718: UserWarning: Data array used to create a new image contains 64-bit ints. This is likely due to creating the array with numpy and passing `int` as the `dtype`. Many tools such as FSL and SPM cannot deal with int64 in Nifti images, so for compatibility the data has been converted to int32. return self.nifti_masker.inverse_transform(self.data) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/plotting.py:135: UserWarning: Casting data from int32 to float32 return view_img( /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/plotting.py:135: UserWarning: Resampling binary images with continuous or linear interpolation. This might lead to unexpected results. You might consider using nearest interpolation instead. return view_img(
To use any mask we just need to index it by the correct label.
Let's start by using the vmPFC mask (ROI=32) to use as a seed in a functional connectivity analysis.
vmpfc = smoothed.extract_roi(mask=mask_x[32])
_fig = go.Figure()
_fig.add_trace(go.Scatter(
y=vmpfc, mode='lines', name='vmPFC',
line=dict(width=2),
hovertemplate='TR=%{x}<br>intensity=%{y:.2f}<extra></extra>',
))
_fig.update_layout(
xaxis_title='Time (TRs)', yaxis_title='Mean Intensity',
height=350, hovermode='x unified',
margin=dict(l=60, r=20, t=20, b=50),
)
_fig
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:718: UserWarning: Data array used to create a new image contains 64-bit ints. This is likely due to creating the array with numpy and passing `int` as the `dtype`. Many tools such as FSL and SPM cannot deal with int64 in Nifti images, so for compatibility the data has been converted to int32. return self.nifti_masker.inverse_transform(self.data)
Okay, now let's build our regression design matrix to perform the whole-brain functional connectivity analysis.
The goal is to find which regions in the brain have a similar time course to the vmPFC, controlling for all of our covariates (i.e., nuisance regressors).
Functional connectivity analyses are particularly sensitive to artifacts that might induce a temporal relationship, particularly head motion (See this article by Jonathan Power for more details). This means that we will need to use slightly different steps to preprocess data for this type of analyis then a typical event related mass univariate analysis.
We are going to remove the mean from our vmPFC signal. We are also going to include the average activity in CSF as an additional nuisance regressor to remove physiological artifacts. Finally, we will be including our 24 motion covariates as well as linear and quadratic trends. We need to be a little careful about filtering as the normal high pass filter for an event related design might be too short and will remove potential signals of interest.
Resting state researchers also often remove the global signal, which can reduce physiological and motion related artifacts and also increase the likelihood of observing negative relationships with your seed regressor (i.e., anticorrelated). This procedure has remained quite controversial in practice (see here here, here, and here for a more in depth discussion). We think that in general including covariates like CSF should be sufficient. It is also common to additionally include covariates from white matter masks, and also multiple principal components of this signal rather than just the mean (see more details about compcorr.
Overall, this code should seem very familiar as it is pretty much the same procedure we used in the single subject GLM tutorial. However, instead of modeling the task design, we are interested in calculating the functional connectivity with the vmPFC.
tr = get_tr()
_fwhm = 6
n_tr = len(data)
def make_motion_covariates(mc, tr):
z_mc = zscore(mc)
parts = {
'': z_mc,
'_sq': z_mc ** 2,
'_diff': z_mc.diff(),
'_diff_sq': z_mc.diff() ** 2,
}
all_mc = pd.concat(
[df.rename(columns=lambda c: f'{c}{suffix}') for suffix, df in parts.items()],
axis=1,
)
all_mc.fillna(value=0, inplace=True)
return Design_Matrix(all_mc, sampling_freq=1 / tr)
vmpfc_1 = zscore(pd.DataFrame(vmpfc, columns=['vmpfc']))
_csf_mask = Brain_Data(get_csf_mask_path(sub)).threshold(upper=0.7, binarize=True)
csf = zscore(pd.DataFrame(smoothed.extract_roi(mask=_csf_mask).T, columns=['csf']))
spikes = smoothed.find_spikes(global_spike_cutoff=3, diff_spike_cutoff=3)
_covariates = load_confounds(sub)
_mc = _covariates[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']]
mc_cov = make_motion_covariates(_mc, tr)
dm = Design_Matrix(pd.concat([vmpfc_1, csf, mc_cov, spikes.drop(labels='TR', axis=1)], axis=1), sampling_freq=1 / tr)
dm = dm.add_poly(order=2, include_lower=True)
dm.convolved = ['vmpfc']
smoothed.X = dm
_stats_vmpfc = smoothed.regress()
vmpfc_conn = _stats_vmpfc['beta']
Warning: You are sending unauthenticated requests to the HF Hub. Please set a HF_TOKEN to enable higher rate limits and faster downloads. /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:253: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = self.nifti_masker.fit_transform(data) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:253: UserWarning: imgs are being resampled to the mask_img resolution. This process is memory intensive. You might want to provide a target_affine that is equal to the affine of the imgs or resample the mask beforehand to save memory and computation time. self.data = self.nifti_masker.fit_transform(data) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/stats.py:1789: ChainedAssignmentError: A value is being set on a copy of a DataFrame or Series through chained assignment. Such chained assignment never works to update the original DataFrame or Series, because the intermediate object on which we are setting values always behaves as a copy (due to Copy-on-Write). Try using '.loc[row_indexer, col_indexer] = value' instead, to perform the assignment in a single step. See the documentation for a more detailed explanation: https://pandas.pydata.org/pandas-docs/stable/user_guide/copy_on_write.html#chained-assignment outlier["diff_spike" + str(i + 1)].iloc[int(loc)] = 1 /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:932: UserWarning: If design matrices are supplied, [drift_model, t_r] will be ignored. glm.fit(self.to_nifti(), design_matrices=self.X) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data))
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order)
Notice how this analysis identifies the default network? This analysis is very similar to the original papers that identified the default mode network using resting state data.
For an actual analysis, we would need to repeat this procedure over all of the participants in our sample and then perform a second level group analysis to identify which voxels are consistently coactive with the vmPFC. We will explore group level analyses in the exercises.
Psychophysiological Interactions¶
Suppose we were interested in seeing if the vmPFC was connected to other regions differently when performing a finger tapping task compared to all other conditions. To compute this analysis, we will need to create a new design matrix that combines the motor regressors and then calculates an interaction term between the seed region activity (e.g., vmpfc) and the condition of interest (e.g., motor).
This type of analysis called, psychophysiological interactions was originally proposed by Friston et al., 1997. For a more hands on and practical discussion read this paper and watch this video by Jeanette Mumford and a follow up video) of a more generalized method.
_img = nib.load(get_file('S01', 'derivatives', 'bold'))
print(f'shape: {_img.shape}')
print(f'voxel size: {_img.header.get_zooms()}')
print(f'TR: {_img.header.get_zooms()[3]} s')
_img
shape: (91, 109, 91, 128) voxel size: (np.float32(2.0), np.float32(2.0), np.float32(2.0), np.float32(1.0)) TR: 1.0 s
<nibabel.nifti1.Nifti1Image object at 0x7f75ee8f8950>
def load_bids_events(subject):
"""Create a Design_Matrix from BIDS event file.
Bypasses nltools.onsets_to_dm because that wrapper has a column-name
mangling bug when TR is supplied. Calls nilearn directly with
hrf_model='glover', then resets the index so the result can be
horizontally concatenated with other 0..n_tr-indexed DataFrames.
"""
from nilearn.glm.first_level import make_first_level_design_matrix as _make_dm
tr = get_tr()
n_tr = nib.load(get_file(subject, 'derivatives', 'bold')).shape[-1]
onsets = load_events(subject)
frame_times = np.arange(n_tr) * tr
dm_raw = _make_dm(frame_times, events=onsets, hrf_model='glover', drift_model=None)
convolved = [c for c in dm_raw.columns if c != 'constant']
return Design_Matrix(dm_raw, convolved=convolved, sampling_freq=1 / tr,
polys=['constant']).reset_index(drop=True)
dm_1 = load_bids_events('S01')
motor_variables = ['video_left_hand', 'audio_left_hand', 'video_right_hand', 'audio_right_hand']
ppi_dm = dm_1.drop(motor_variables, axis=1)
ppi_dm['motor'] = pd.Series(dm_1.loc[:, motor_variables].sum(axis=1))
ppi_dm_conv = ppi_dm.convolve()
ppi_dm_conv['vmpfc'] = vmpfc_1.values
ppi_dm_conv['vmpfc_motor'] = ppi_dm_conv['vmpfc'] * ppi_dm_conv['motor_c0']
dm_1 = Design_Matrix(pd.concat([ppi_dm_conv, csf, mc_cov, spikes.drop(labels='TR', axis=1)], axis=1), sampling_freq=1 / tr)
dm_1 = dm_1.add_poly(order=2, include_lower=True)
dm_1.convolved = list(ppi_dm_conv.columns)
dm_1.heatmap()
plt.clf()
Okay, now we are ready to run the regression analysis and inspect the interaction term to find regions where the connectivity profile changes as a function of the motor task.
We will run the regression and smooth all of the images, and then examine the beta image for the PPI interaction term.
smoothed.X = dm_1
ppi_stats = smoothed.regress()
vmpfc_motor_ppi = ppi_stats['beta'][int(np.where(smoothed.X.columns == 'vmpfc_motor')[0][0])]
vmpfc_motor_ppi.iplot()
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:932: UserWarning: If design matrices are supplied, [drift_model, t_r] will be ignored. glm.fit(self.to_nifti(), design_matrices=self.X) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data))
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/plotting.py:135: UserWarning: The given float value must not exceed 3.84581147309937e-13. But, you have given threshold=1e-06. return view_img( /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/plotting.py:135: UserWarning: The given float value must not exceed 0.0. But, you have given threshold=1e-06. return view_img( /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order)
This analysis tells us which regions are more functionally connected with the vmPFC during the motor conditions relative to the rest of experiment.
We can make a thresholded interactive plot to interrogate these results, but it looks like it identifies the ACC/pre-SMA.
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/plotting.py:135: UserWarning: The given float value must not exceed 3.84581147309937e-13. But, you have given threshold=1e-06. return view_img( /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/plotting.py:135: UserWarning: The given float value must not exceed 0.0. But, you have given threshold=1e-06. return view_img( /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order)
Dynamic Connectivity¶
All of the methods we have discussed so far assume that the relationship between two regions is stationary - or remains constant over the entire dataset. However, it is possible that voxels are connected to other voxels at specific points in time, but then change how they are connected when they are computing a different function or in different psychological state.
Time-varying connectivity is beyond the scope of the current tutorial, but we encourage you to watch this video from Principles of fMRI for more details
Effective Connectivity¶
Effective connectivity refers to the degree that one brain region has a directed influence on another region. This approach requires making a number of assumptions about the data and requires testing how well a particular model describes the data. Typically, most researchers will create a model of a small number of nodes and compare different models to each other. This is because the overall model fit is typically in itself uninterpretable and because formulating large models can be quite difficult and computationally expensive. The number of connections can be calculated as:
\(connections = \frac{n(n-1)}{2}\), where \(n\) is the total number of nodes.
Let's watch a short video by Martin Lindquist that provides an overview to different approaches to effectivity connectivity.
Structural Equation Modeling
Structural equation modeling (SEM) is one early technique that was used to model the causal relationship between multiple nodes. SEM requires specifying a causal relationship between nodes in terms of a set of linear equations. The parameters of this system of equations reflects the connectivity matrix. Users are expected to formulate their own hypothesized relationship between variables with a value of one when there is an expected relationship, and zero when there is no relationship. Then we estimate the parameters of the model and evaluate how well the model describes the observed data.
We will not be discussing this method in much detail. In practice, this method is more routinely used to examine how brain activations mediate relationships between other regions, or between different psychological constructs (e.g., X -> Z -> Y).
Here are a couple of videos specifically examining how to conduct mediation and moderation analyses from Principles of fMRI (Mediation and Moderation Part I,
Mediation and Moderation Part II)Granger Causality¶
Granger causality was originally developed in econometrics and is used to determine temporal causality. The idea is to quantify how past values of one brain region predict the current value of another brain region. This analysis can also be performed in the frequency domain using measures of coherence between two regions. In general, this technique is rarely used in fMRI data analysis as it requires making assumptions that all regions have the same hemodynamic response function (which does not seem to be true), and that the relationship is stationary, or not varying over time.
Here is a video from Principles of fMRI explaining Granger Causality in more detail.
Dynamic Causal Modeling¶
Dynamic Causal Modeling (DCM) is a method specifically developed for conducting causal analyses between regions of the brain for fMRI data. The key innovation is that the developers of this method have specified a generative model for how neuronal firing will be reflected in observed BOLD activity. This addresess one of the problems with SEM, which assumes that each ROI has the same hemodynamic response.
In practice, DCM is computationally expensive to estimate and researchers typically specify a couple small models and perform a model comparison (e.g., bayesian model comparison) to determine, which model best explains the data from a set of proposed models.
Here is a video from Principles of fMRI explaining Dynamic Causal Modeling in more detail.
(multivariate-decomposition)=
Multivariate Decomposition¶
So far we have discussed functional connectivity in terms of pairs of regions. However, voxels are most likely not independent from each other and we may want to figure out some latent spatial components that are all functionally connected with each (i.e., covary similarly in time).
To do this type of analysis, we typically use what are called multivariate decomposition methods, which attempt to factorize a data set (i.e., time by voxels) into a lower dimensional set of components, where each has their own unique time course.
The most common decomposition methods or Principal Components Analysis (PCA) and Independent Components Analysis (ICA). You may recall that we have already played with ICA in the beginning of the course in the ICA lab.
Let's learn more about the details of decomposition by watching a short video by Martin Lindquist.
Principal Components Analysis
Principal Components Analysis (PCA) is a multivariate procedure that attempts to explain the variance-covariance structure of a high dimensional random vector. In this procedure, a set of correlated variables are transformed int a set of uncorrelated variables, ordered by the amount of variance in the data that they explain. In fMRI, we use PCA to find spatial maps or eigenimages in the data. This is usually computed using Singular Value Decomposition (SVD). This operation is defined as:
SVD is conceptually very similar to regression. We are trying to explain a matrix _csf_mask = Brain_Data(get_csf_mask_path(sub)).threshold(upper=0.7, binarize=True)
csf_1 = zscore(pd.DataFrame(smoothed.extract_roi(mask=_csf_mask).T, columns=['csf']))
spikes_1 = smoothed.find_spikes(global_spike_cutoff=3, diff_spike_cutoff=3)
_covariates = load_confounds(sub)
_mc = _covariates[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']]
mc_cov_1 = make_motion_covariates(_mc, tr)
dm_2 = Design_Matrix(pd.concat([vmpfc_1, csf_1, mc_cov_1, spikes_1.drop(labels='TR', axis=1)], axis=1), sampling_freq=1 / tr)
dm_2 = dm_2.add_poly(order=2, include_lower=True)
dm_2.convolved = ['vmpfc']
smoothed.X = dm_2
_stats_denoise = smoothed.regress()
smoothed_denoised = _stats_denoise['residual']
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:253: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = self.nifti_masker.fit_transform(data) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:253: UserWarning: imgs are being resampled to the mask_img resolution. This process is memory intensive. You might want to provide a target_affine that is equal to the affine of the imgs or resample the mask beforehand to save memory and computation time. self.data = self.nifti_masker.fit_transform(data) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/stats.py:1789: ChainedAssignmentError: A value is being set on a copy of a DataFrame or Series through chained assignment. Such chained assignment never works to update the original DataFrame or Series, because the intermediate object on which we are setting values always behaves as a copy (due to Copy-on-Write). Try using '.loc[row_indexer, col_indexer] = value' instead, to perform the assignment in a single step. See the documentation for a more detailed explanation: https://pandas.pydata.org/pandas-docs/stable/user_guide/copy_on_write.html#chained-assignment outlier["diff_spike" + str(i + 1)].iloc[int(loc)] = 1 /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:932: UserWarning: If design matrices are supplied, [drift_model, t_r] will be ignored. glm.fit(self.to_nifti(), design_matrices=self.X) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data)) /home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:262: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used. self.data = np.array(self.nifti_masker.fit_transform(data))
Now let's run a PCA on this participant's denoised data. To do this, we will use the .decompose() method from nltools. All we need to do is specify the algorithm we want to use, the dimension we want to reduce (i.e., time - 'images' or space 'voxels'), and the number of components to estimate. Usually, we will be looking at reducing space based on similarity in time, so we will set axis='images'.
n_components = 10
pca_stats_output = smoothed_denoised.decompose(algorithm='pca', axis='images', n_components=n_components)
Now let's inspect the components with our interactive component viewer. Remember the ICA tutorial? Hopefully, you are now better able to understand everything.
component_slider = mo.ui.slider(
start=0, stop=len(pca_stats_output['components']) - 1, value=0, step=1,
label='Component', show_value=True, full_width=True,
)
threshold_slider = mo.ui.slider(
start=0.0, stop=4.0, value=2.0, step=0.1,
label='Threshold (z)', show_value=True, full_width=True,
)
mo.hstack([component_slider, threshold_slider], justify='start')
_comp = component_slider.value
_thresh = threshold_slider.value
# z-score and threshold the selected component map
_component = pca_stats_output['components'][_comp]
_zscored = (_component - _component.mean()) * (1 / _component.std())
_zscored.data = np.where(np.abs(_zscored.data) <= _thresh, 0, _zscored.data)
# Title with variance explained (if PCA)
import sklearn.decomposition as _skd
_decomp = pca_stats_output['decomposition_object']
if isinstance(_decomp, _skd.PCA):
_var = _decomp.explained_variance_ratio_[_comp]
_title = f'Component {_comp}/{len(pca_stats_output["components"])} — Variance Explained: {_var:.2%}'
else:
_title = f'Component {_comp}/{len(pca_stats_output["components"])}'
# Timecourse + power spectrum (plotly subplots)
_weights = pca_stats_output['weights'][:, _comp]
_y = fft(_weights)
_freqs = fftfreq(len(_y), d=get_tr())
_pos_mask = _freqs > 0
_fig = make_subplots(
rows=2, cols=1, vertical_spacing=0.18,
subplot_titles=(f'Timecourse (TR={get_tr()}s)', 'Power Spectrum'),
)
_fig.add_trace(go.Scatter(
y=_weights, mode='lines', line=dict(color='red', width=2),
hovertemplate='TR=%{x}<br>intensity=%{y:.3f}<extra></extra>',
name='timecourse',
), row=1, col=1)
_fig.add_trace(go.Scatter(
x=_freqs[_pos_mask], y=np.abs(_y)[_pos_mask] ** 2,
mode='lines', line=dict(width=2),
hovertemplate='f=%{x:.3f} Hz<br>power=%{y:.2f}<extra></extra>',
name='power',
), row=2, col=1)
_fig.update_xaxes(title_text='Time (TRs)', row=1, col=1)
_fig.update_yaxes(title_text='Intensity (AU)', row=1, col=1)
_fig.update_xaxes(title_text='Frequency (Hz)', row=2, col=1)
_fig.update_yaxes(title_text='Power', row=2, col=1)
_fig.update_layout(
title=_title, height=600, showlegend=False,
margin=dict(l=60, r=20, t=80, b=50),
)
mo.vstack([_zscored.iplot(), _fig])
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order)
We can also examine the eigenvalues/singular values or scaling factor of each, which are the diagonals of \(S\).
These values are stored in the 'decomposition_object' of the stats_output and are in the variable called .singular_values_.
_sv = pca_stats_output['decomposition_object'].singular_values_
_fig = go.Figure()
_fig.add_trace(go.Scatter(
x=np.arange(len(_sv)), y=_sv, mode='lines+markers',
line=dict(width=2),
hovertemplate='component %{x}<br>singular value=%{y:.3f}<extra></extra>',
))
_fig.update_layout(
xaxis_title='Component', yaxis_title='Singular Value',
height=350, margin=dict(l=60, r=20, t=20, b=50),
)
_fig
We can use these values to calculate the overall variance explained by each component. These values are stored in the 'decomposition_object' of the stats_output and are in the variable called .explained_variance_ratio_.
These values can be used to create what is called a scree plot to figure out the percent variance of \(X\) explained by each component. Remember, in PCA, components are ordered by descending variance explained.
_var = pca_stats_output['decomposition_object'].explained_variance_ratio_
_cum = np.cumsum(_var)
_fig = make_subplots(rows=1, cols=2, subplot_titles=('Variance Explained', 'Cumulative Variance Explained'))
_fig.add_trace(go.Scatter(x=np.arange(len(_var)), y=_var, mode='lines+markers',
name='per-component', line=dict(width=2),
hovertemplate='component %{x}<br>variance=%{y:.3%}<extra></extra>'),
row=1, col=1)
_fig.add_trace(go.Scatter(x=np.arange(len(_cum)), y=_cum, mode='lines+markers',
name='cumulative', line=dict(width=2),
hovertemplate='component %{x}<br>cumulative=%{y:.3%}<extra></extra>'),
row=1, col=2)
_fig.update_xaxes(title_text='Component', row=1, col=1)
_fig.update_xaxes(title_text='Component', row=1, col=2)
_fig.update_yaxes(title_text='Percent Variance Explained', row=1, col=1, tickformat='.0%')
_fig.update_yaxes(title_text='Percent Variance Explained', row=1, col=2, tickformat='.0%')
_fig.update_layout(height=400, showlegend=False, margin=dict(l=60, r=20, t=50, b=50))
_fig
Independent Components Analysis¶
Independent Components Analysis (ICA) is a method to blindly separate a source signal into spatially independent components. This approach assumes that the data consts of \(p\) spatially independent components, which are linearly mixed and spatially fixed. PCA assumes orthonormality constraint, while ICA only assumes independence.
\(X = AS\), where \(A\) is the mixing matrix and \(S\) is the source matrix
In ICA we find an un-mixing matrix \(W\), such that \(Y = WX\) provides an approximation to \(S\). To estimate the mixing matrix, ICA assumes that the sources are (1) linearly mixed, (2) the components are statistically independent, and (3) the components are non-Gaussian.
It is trivial to run ICA on our data as it only requires switching algorithm='pca' to algorithm='ica' when using the decompose() method.
We will experiment with this in our exercises.
We also encourage interested readers to check out our more in depth ICA tutorial.
Graph Theory
Similar to describing the structure of social networks, graph theory has also been used to characterize regions of the brain based on how they connected to other regions. Nodes in the network typically describe specific brain regions and edges represent the strength of the association between each edge. That is, the network can be represented as a graph of pairwise relationships between each region of the brain. There are many different metrics of graphs that can be used to describe the overall efficiency of a network (e.g., small worldness), or how connected a region is to other regions (e.g., degree, centrality), or how long it would take to send information from one node to another node (e.g., path length, connectivity).
Let's watch a short video by Martin Lindquist providing a more in depth introduction to graph theory.Suppose, we were interested in identifying which regions of the brain had the highest degree of centrality based on functional connectivity. There are many different ways to do this, but they all involve specifying a set of nodes (i.e., ROIs) and calculating the edges between each node. Finally, we would need to pick a centrality metric and calculate the overall level of centrality for each region.
Let's do this quickly building off of our seed-based functional connectivity analysis.
Similar, to the PCA example, let's work with the denoised data. First, let's extract the average time course within each ROI from our 50 parcels and plot the results.
rois = smoothed_denoised.extract_roi(mask=mask_1)
_fig = go.Figure()
for i in range(rois.shape[0]):
_fig.add_trace(go.Scatter(
y=rois[i], mode='lines', name=f'ROI {i}', line=dict(width=1),
hovertemplate=f'ROI {i}<br>TR=%{{x}}<br>intensity=%{{y:.2f}}<extra></extra>',
))
_fig.update_layout(
xaxis_title='Time (TRs)', yaxis_title='Mean Intensity',
height=400, showlegend=False, hovermode='closest',
margin=dict(l=60, r=20, t=20, b=50),
)
_fig
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/nltools/data/brain_data.py:718: UserWarning: Data array used to create a new image contains 64-bit ints. This is likely due to creating the array with numpy and passing `int` as the `dtype`. Many tools such as FSL and SPM cannot deal with int64 in Nifti images, so for compatibility the data has been converted to int32. return self.nifti_masker.inverse_transform(self.data)
Now that we have specified our 50 nodes, we need to calculate the edges of the graph. We will be using pearson correlations. We will be using the pairwise_distances function from scikit-learn as it is much faster than most other correlation measures. We will then convert the distance metric into similarities by subtracting all of the values from 1.
Let's visualize the resulting correlation matrix as a heatmap using seaborn.
roi_corr = 1 - pairwise_distances(rois, metric='correlation')
_fig = go.Figure(go.Heatmap(
z=roi_corr, zmin=-1, zmax=1, colorscale='RdBu_r', reversescale=True,
hovertemplate='ROI %{x} → %{y}<br>r=%{z:.3f}<extra></extra>',
))
_fig.update_layout(
title='ROI-to-ROI correlation', height=600, width=600,
yaxis=dict(autorange='reversed', scaleanchor='x'),
margin=dict(l=60, r=20, t=50, b=50),
)
_fig
Now we need to convert this correlation matrix into a graph and calculate a centrality measure. We will use the Adjacency class from nltools as it has many functions that are useful for working with this type of data, including casting these type of matrices into networkx graph objects.
We will be using the networkx python toolbox to work with graphs and compute different metrics of the graph.
Let's calculate degree centrality, which is the total number of nodes each node is connected with. Unfortunately, many graph theory metrics require working with adjacency matrices, which are binary matrices indicating the presence of an edge or not. To create this, we will simply apply an arbitrary threshold to our correlation matrix.
_a = Adjacency(roi_corr, matrix_type='similarity', labels=[x for x in range(50)])
a_thresholded = _a.threshold(upper=0.6, binarize=True)
a_thresholded.plot()
Okay, now that we have a thresholded binary matrix, let's cast our data into a networkx object and calculate the degree centrality of each ROI and make a quick plot of the graph.
G = a_thresholded.to_graph()
_pos = nx.kamada_kawai_layout(G)
_node_degree = dict(G.degree())
# Edges: flatten to (x0, x1, None, x0, x1, None, ...) for a single Scatter trace
_edge_x = []
_edge_y = []
for u, v in G.edges():
x0, y0 = _pos[u]
x1, y1 = _pos[v]
_edge_x.extend([x0, x1, None])
_edge_y.extend([y0, y1, None])
# Nodes
_node_x = [_pos[n][0] for n in G.nodes()]
_node_y = [_pos[n][1] for n in G.nodes()]
_node_text = [f'ROI {n}<br>degree={_node_degree[n]}' for n in G.nodes()]
_node_labels = [str(n) for n in G.nodes()]
_node_sizes = [max(_node_degree[n] * 4, 8) for n in G.nodes()]
_node_colors = [_node_degree[n] for n in G.nodes()]
_fig = go.Figure()
_fig.add_trace(go.Scatter(
x=_edge_x, y=_edge_y, mode='lines',
line=dict(width=1.5, color='rgba(80,80,100,0.25)'),
hoverinfo='skip', showlegend=False,
))
_fig.add_trace(go.Scatter(
x=_node_x, y=_node_y, mode='markers+text',
marker=dict(
size=_node_sizes, color=_node_colors, colorscale='Reds',
reversescale=True, line=dict(color='darkslategray', width=2),
colorbar=dict(title='Degree', thickness=15),
),
text=_node_labels, textposition='middle center',
textfont=dict(color='darkslategray', size=10),
hovertext=_node_text, hoverinfo='text', showlegend=False,
))
_fig.update_layout(
height=700, width=900,
xaxis=dict(showgrid=False, zeroline=False, showticklabels=False),
yaxis=dict(showgrid=False, zeroline=False, showticklabels=False, scaleanchor='x'),
margin=dict(l=20, r=20, t=20, b=20),
plot_bgcolor='white',
)
_fig
We can also plot the distribution of degree using this threshold.
_degrees = list(dict(G.degree).values())
_fig = go.Figure(go.Histogram(
x=_degrees, name='degree',
hovertemplate='degree=%{x}<br>count=%{y}<extra></extra>',
))
_fig.update_layout(
xaxis_title='Degree', yaxis_title='Frequency',
height=350, margin=dict(l=60, r=20, t=20, b=50),
)
_fig
What if we wanted to map the degree of each node back onto the brain?
This would allow us to visualize which of the parcels had more direct pairwise connections.
To do this, we will simply scale our expanded binary mask object by the node degree. We will then combine the masks by concatenating through recasting as a brain_data object and then summing across all ROIs.
degree = pd.Series(dict(G.degree()))
brain_degree = roi_to_brain(degree, mask_x)
brain_degree.iplot()
[2K [2K
/home/runner/work/dartbrains/dartbrains/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:840: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray. a.partition(kth, axis=axis, kind=kind, order=order)
This analysis shows that the insula is one of the regions that appears to have the highest degree in this analysis. This is a fairly classic finding with the insula frequently found to be highly connected with other regions. Of course, we are only looking at one subject in a very short task (and selecting a completely arbitrary cutoff). We would need to show this survives correction after performing a group analysis.
Exercises¶
Let's practice what we learned through a few different exercises.
1) Let's calculate seed-based functional connectivity using a different ROI - the right motor cortex¶
- Calculate functional connectivity using roi=48 with the whole brain.
2) Calculate a group level analysis for this connectivity analysis¶
- this will require running this analysis over all subjects
- then running a one sample t-test
- then correcting for multiple tests with fdr.
3) Calculate an ICA¶
- run an ICA analysis for subject01 with 5 components
- plot each spatial component and its associated timecourse
4) Calculate Eigenvector Centrality for each Region¶
- figure out how to calculate eigenvector centrality and compute it for each region.
5) Calculate a group level analysis for this graph theoretic analysis¶
- this will require running this analysis over all subjects
- then running a one sample t-test
- then correcting for multiple tests with fdr.