State-based dFC for voxelwise analysis

Voxelwise state-based dFC analysis is supported naturally by reshaping data to 2D (unit, time) and feeding it into the same algorithms used for ROI-level analyses.

• Units are voxels (or surface vertices).

• Time is the number of volumes / TRs.

For example, with a NIfTI and brain mask, we can simply transform it into a 2D array:

import numpy as np
from nilearn import datasets, maskers, plotting, image
from comet import connectivity, utils

# This will take 1-2 minutes to download
save_dir = "."
haxby = datasets.fetch_haxby(data_dir=save_dir, subjects=1)

# Get fMRI and mask paths
fmri_path = haxby.func[0]
mask_path = haxby.mask_vt[0] # Chose the ventral-temporal mask

# Extract voxelwise time series array
masker = maskers.NiftiMasker(mask_img=mask_path, standardize="zscore_sample", detrend=True, dtype='float32').fit()
time_series = masker.transform(fmri_path)
print("Time series shape:", time_series.shape)

[fetch_haxby] Dataset found in haxby2001

/tmp/ipykernel_95734/25630856.py:14: UserWarning: The provided image has no sform in its header. Please check the provided file. Results may not be as expected.
  masker = maskers.NiftiMasker(mask_img=mask_path, standardize="zscore_sample", detrend=True, dtype='float32').fit()

Time series shape: (1452, 577)

We can then estimate the coactivation patterns:

cap = connectivity.CoactivationPatterns(time_series=time_series, n_states=5)
state_tc, states = cap.estimate()

print("state_tc shape:", state_tc.shape)
print("states shape:  ", states.shape)

fig1, ax1 = utils.state_plots(states=states, figsize=(8,2))
fig2, ax2 = utils.state_plots(state_tc=state_tc, figsize=(8,2))

CAP: 100%|██████████| 1/1 [00:00<00:00,  1.37it/s]

state_tc shape: (1, 1452)
states shape:   (577, 577, 5)

The states are (\(voxels * voxels\)) connectivity matrices created as outer products of the group centroids. To obtain a voxelwise activation pattern per state for visualisation, we extract the top eigenvector of each state matrix (positive orientation):

def top_eigenvector(A):
    # Compute leading eigenpair
    vals, vecs = np.linalg.eigh(A)
    v = vecs[:, -1]

    # Make the largest-magnitude element positive
    if np.abs(v).max() > 0 and np.sign(v[np.argmax(np.abs(v))]) < 0:
        v = -v
    return v.astype(np.float32)

state_vectors = np.stack([top_eigenvector(states[:, :, s]) for s in range(states.shape[2])], axis=0)  # (n_states, P)
print("Recovered state vectors (n_states, P):", state_vectors.shape)

Recovered state vectors (n_states, P): (5, 577)

We can invert the vectors back to 3D images to plot the coactivation pattern:

img = masker.inverse_transform(state_vectors)  # shape: (x, y, z, n_states)
plotting.plot_stat_map(image.index_img(img, 1), display_mode="ortho", title=f"CAP State 1", cmap="coolwarm");