Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.

Code ↔ Paper

18 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 18 matches

[1] § Methods › Simulation of phase delay structure in BOLD-like 1/f signals ↔ BOLD_simulation.py, lines 100–160 · score 0.82 · Wilson Cowan, 1–100 s, OU, field, simulated, 1 s
[2] § Results › Behavioral phenotypes are associated with CPCA-derived features ↔ sCCA_analysis_wholebrain.py, lines 349–498 · score 0.82 · aggressive behavior, language comprehension, working memory, canonical components, fear, FDR
[3] § Methods › Comparison with traditional spatiotemporal models ↔ figures.ipynb, lines 1580–1725 · score 0.78 · temporal ICA, spatial ICA, BOLD activity, spatial weight, spatial correlations, TICA
[4] § Methods › Comparison with traditional spatiotemporal models › Quasi-periodic patterns (QPPs) ↔ run_qpp.py, lines 30–127 · score 0.77 · window length, correlation threshold, QPPs, iteratively, segment, scan
[5] § Methods › Behavior data ↔ sCCA_analysis_wholebrain.py, lines 349–498 · score 0.77 · aggressive behavior, language comprehension, working memory, fear, alertness, personality
[6] § Methods › Sex classification ↔ sex_LSVM.py, lines 324–392 · score 0.77 · hidden layer, confusion matrices, ROC, fold, training, curves
[7] § Results › CPCA-based reconstruction preserves functional connectivity modularity ↔ figures.ipynb, lines 1203–1293 · score 0.76 · Louvain community detection, module partition, original correlation, NMI, functional connectivity, reconstructed
[8] § Methods › Signal reconstruction and functional connectivity ↔ figures.ipynb, lines 1203–1293 · score 0.71 · normalized mutual information, original FC, Louvain, NMI, reconstructed, matrices
[9] § Results › Three prominent spatiotemporal patterns of the cerebellum reflect FC topographies ↔ run_analysis.sh, lines 1–63 · score 0.68 · hidden Markov models, temporal ICA, spatial ICA, HMM, eigenmaps, FC
[10] § Methods › Comparison with traditional spatiotemporal models ↔ CPCA_analysis.py, lines 400–467 · score 0.66 · voxel coordinates, MNI, Moran, Pearson, HMM, templates
[11] § Methods › Statistical analysis ↔ CPCA_analysis.py, lines 182–203 · score 0.66 · Moran spectral randomization, spatial autocorrelation, Permutation, connectivity, components
[12] § Results › Sex differences are detected by machine learning methods ↔ ROC_CM_null.py, lines 124–158 · score 0.64 · ROC curve, confusion matrix, sex classification, ANN, accuracy, predictions
[13] § Results › Three prominent spatiotemporal patterns of the cerebellum reflect FC topographies ↔ figures.ipynb, lines 724–766 · score 0.60 · hidden Markov models, FC topographies, TICA, SICA, zero lag, HMM
[14] § Results › Sex differences are detected by machine learning methods ↔ sex_LSVM.py, lines 324–392 · score 0.59 · ROC curve, confusion matrix, AUC, ANN, LSVM, accuracy
[15] § Methods › Sparse canonical correlation analysis ↔ sCCA_analysis_wholebrain.py, lines 37–185 · score 0.57 · sparse canonical correlation, sCCA, Model
[16] § Methods › Sex classification ↔ ROC_CM_null.py, lines 124–158 · score 0.56 · confusion matrices, ROC, curves, ANN, threshold, accuracy
[17] § Results › Phase delay components capture temporal propagation patterns ↔ figures.ipynb, lines 1580–1725 · score 0.55 · identified recurring, low dimensional, BOLD activity, 0.01 Hz, temporally, propagation
[18] § Results ↔ sCCA_analysis_wholebrain.py, lines 545–596 · score 0.51 · sCCA, rsFC, subsets, sparse, trained, transform

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 2,209 lines · 100 KB · no license · 5 matches

# %% [markdown]
# # Study Abstract
#
# Low-frequency blood-oxygenated-level-dependent (BOLD) signals are known to exhibit a global spatiotemporal pattern of activity at rest, labeled the quasiperiodic pattern (QPP). This spatiotemporal pattern can be described as a recurring global propagation of BOLD activity from task-positive network regions towards default mode network regions. Alongside the QPP, a wide variety of functional connectivity topographies contrasting task-positive and default mode network regions have been reported in the literature. In this study, we demonstrate that most widely-studied functional connectivity topographies arise from the same spatiotemporal dynamics - the QPP. In other words, much of the resting-state fMRI literature has been describing the same thing with different methods. Using resting-state functional magnetic resonance imaging scans from the Human Connectome Project (n=50), we examined the relationship between previously observed functional connectivity topographies and the time-lagged dynamics of the QPP. We find that the time-lagged dynamics of the QPP are distributed across multiple axes in low-dimensional functional connectivity space(s). Popular functional connectivity topographies generally correspond to one or more of these axes. Thus, functional connectivity topographies represent low-dimensional descriptions of the synchronous dynamics within the larger time-lagged QPP pattern. We further demonstrate that ‘global signal’ and ‘task-positive/task-negative’ functional connectivity topographies arise from the same dynamics of the QPP: a time-lag in BOLD signals between task-positive and task-negative network regions. Overall, we find that the QPP underlies a striking variety of previously observed phenomena in low-frequency resting-state BOLD signals.
# %% [markdown]
# ## Module Imports
# %%
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import nibabel as nb
import numpy as np
import pandas as pd
import pickle
from brainspace.null_models import SpinPermutations
from bct import community_louvain, partition_distance
from IPython.core.display import HTML
from matplotlib.offsetbox import OffsetImage, AnnotationBbox, TextArea
from matplotlib import cm, colors, animation
from matplotlib.cm import ScalarMappable
from matplotlib.colors import LinearSegmentedColormap, Normalize
from matplotlib.patches import FancyArrowPatch, FancyBboxPatch, BoxStyle, Rectangle
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.mplot3d.axes3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from mpl_toolkits.axes_grid1 import AxesGrid
from nilearn.surface import load_surf_mesh
from numpy import random as rand
from run_main_pca import pca as cpca, rotation
from scipy import linalg
from scipy.signal import resample, welch, hilbert
from scipy.spatial.distance import cdist, pdist, squareform
from scipy.stats import zscore
from sklearn.linear_model import LinearRegression
from sklearn.cluster import KMeans
from sklearn.decomposition import KernelPCA, PCA, FactorAnalysis, FastICA
from sklearn.manifold import SpectralEmbedding
from utils.utils import load_data_and_stack, load_gifti, pull_gifti_data, load_cifti, pull_cifti_data, \
write_to_gifti, write_to_cifti
from utils.rotation import varimax
# %% [markdown]
# ## Helper Functions
# %%
# Global variables
tr = 0.72 # HCP sampling rate
n_vertices = 4801 # Number of vertices in functional scan (after pre-processing)
n_vert_L=2562 # Number of vertices in left cortex (n_vert_R = n_vertices - n_vert_L)
n_ts=60000 # Number of time-points in group-concatenated time courses
# Helper Functions
def axis3d_scaler(x, y, z):
# https://stackoverflow.plex_/questions/30223161/matplotlib-mplot3d-how-to-increase-the-size-of-an-axis-stretch-in-a-3d-plot
"""
stretch axis of 3-dimensional plot
"""
scale=np.diag([x, y, z, 1.0])
scale=scale*(1.0/scale.max())
scale[3,3]=1.0
return scale
def circular_corr(alpha1, alpha2, nanrobust, axis=None):
# Taken from:
# https://github.com/jhamrick/python-snippets/blob/master/snippets/circstats.py
"""
Calculate circulation correlation between two phase time series
"""
if axis is not None and alpha1.shape[axis] != alpha2.shape[axis]:
raise(ValueError, "shape mismatch")
# compute mean directions
if axis is None:
n = alpha1.size
else:
n = alpha1.shape[axis]
#################################################################
c1 = np.cos(alpha1)
c1_2 = np.cos(2*alpha1)
c2 = np.cos(alpha2)
c2_2 = np.cos(2*alpha2)
s1 = np.sin(alpha1)
s1_2 = np.sin(2*alpha1)
s2 = np.sin(alpha2)
s2_2 = np.sin(2*alpha2)
if nanrobust:
sumfunc = lambda x: np.nansum(x, axis=axis)
else:
sumfunc = lambda x: np.sum(x, axis=axis)
num = 4 * (sumfunc(c1*c2) * sumfunc(s1*s2) -
sumfunc(c1*s2) * sumfunc(s1*c2))
den = np.sqrt((n**2 - sumfunc(c1_2)**2 - sumfunc(s1_2)**2) *
(n**2 - sumfunc(c2_2)**2 - sumfunc(s2_2)**2))
rho = num / den
return rho
def complex_motion(g1, g2, g2_phase_shift, y, w, n_samples, n_cycles, n_decay_cycles, var,
decay_phase_shift=0):
"""
Generate (damped) 2-dimensional complex motion between two spatial patterns (g1 and g2). Properties of motion
are controlled by parameters described below.
Parameters
----------
g1: 2-d numpy array
spatial pattern 1
g2: 2-d numpy array
spatial pattern 2
g2_phase_shift: float
phase shift, in radians, to be applied to g2 pattern
y : float
Decay amplitude of complex motion. Decay is sinusoidal, so the complex
motion waxes and wanes within a cycle.
w: float
angular frequency of sine oscillation. Controls the frequency or speed of the oscillation
n_samples: int
number of time samples to generate
n_cycles: int
number of full oscillations to cycle through with specified number of time samples (n_samples)
n_decay_cycles: int
number of decay cycles to cycle through with specified number of time samples (n_samples)
var: float
variance of gaussian noise added to each time point
Returns:
list: complex motion spatial pattern at each time point (n_samples)
"""
g1_vec = g1.flatten()
g2_vec = g2.flatten()
x_len, y_len = g1.shape
t_vec = np.linspace(1,(2*np.pi)*n_cycles, n_samples)
t_decay_vec = 0.5*np.cos(
np.linspace(1,(2*np.pi)*n_decay_cycles, n_samples) + decay_phase_shift
)+0.5
g_anim = []
for t, t_d in zip(t_vec, t_decay_vec):
g_grid = []
for g1, g2 in zip(g1_vec, g2_vec):
g_grid.append(2*np.exp(y*t_d)*(np.sin(w*t)*g1 - np.sin(w*t + g2_phase_shift)*g2))
grid = np.array(g_grid).reshape(x_len, y_len)
grid += np.random.normal(0,var,(grid.shape[0], grid.shape[1]))
g_anim.append(grid)
return g_anim
def complex_motion_animation(g_anim, fig, ax, vmin=-1, vmax=1):
"""Play animation of complex motion"""
ims = []
for g in g_anim:
im = ax.imshow(g, cmap='coolwarm', vmin=vmin, vmax=vmax)
ims.append([im])
anim = animation.ArtistAnimation(fig, ims, interval=100, blit=True, repeat_delay=1000)
return anim
def convert_polar_xticks_to_radians_secs(ax, cycle_length):
# https://stackoverflow.com/questions/21172228/python-matplotlib-polar-plots-with-angular-labels-in-radians
# Converts x-tick labels from degrees to radians
# Get the x-tick positions (returns in radians)
label_positions = ax.get_xticks()
# Convert to a list since we want to change the type of the elements
labels = list(label_positions)
# Format each label (edit this function however you'd like)
labels = [format_secs_radians_label(label, cycle_length) for label in labels]
ax.set_xticklabels(labels)
def convert_polar_xticks_to_radians(ax):
# https://stackoverflow.com/questions/21172228/python-matplotlib-polar-plots-with-angular-labels-in-radians
"""
Converts x-tick labels from degrees to radians
"""
# Get the x-tick positions (returns in radians)
label_positions = ax.get_xticks()
# Convert to a list since we want to change the type of the elements
labels = list(label_positions)
# Format each label (edit this function however you'd like)
labels = [format_radians_label(label) for label in labels]
ax.set_xticklabels(labels)
def create_hrf_group(n_ts, activation_indx, ts_len, tr, amplitude, phase_jitter,
amplitude_jitter, ts_sampling=0.01, repeat_n=1):
"""
n_ts: number of timeseries
activation_indx = index of activation time point
ts_len: length of time series
tr: the sampling rate of the original time series
amplitude: amplitude of double gamma function
phase_offset_window: allowable phase offsets between time series -
set as a symmetric window length - sampled from uniform distribution
ts_sampling: resolution of original time series - default=0.01 Hz
std_noise: amount of gaussian noise to add to time series - scaling parameter between 0 and 1
"""
hrf=double_gamma_hrf(60, ts_sampling)
ts_all = np.zeros((n_ts, ts_len))
for n in range(n_ts):
ts = ts_all[n,:]
indx = rand.randint(activation_indx - phase_jitter,
activation_indx + phase_jitter)
amp = rand.randint(amplitude - amplitude_jitter,
amplitude + amplitude_jitter)
ts[indx] = 1
ts_all[n,:] = (convolve_hrf_events(hrf, ts) * amplitude)
n_resample=np.int(ts_sampling*ts_len/tr)
hrf_ts_resample = resample(ts_all, n_resample, axis=1)
return np.tile(hrf_ts_resample, repeat_n)
def cropImage(img, width_l, width_r, height):
"""
Crop image by width (tuple) and height (tuple) ranges
"""
slice1 = img[height[0]:height[1],width_l[0]:width_l[1],:]
slice2 = img[height[0]:height[1],width_r[0]:width_r[1],:]
merged_image = np.append(slice1, slice2, axis=1)
return merged_image
def cropImage_single(img, width, height):
return img[height[0]:-height[1]:,width[0]:-width[1],:]
def format_radians_label(float_in):
# https://stackoverflow.com/questions/21172228/python-matplotlib-polar-plots-with-angular-labels-in-radians
"""
Converts a float value in radians into a string representation of that float
"""
string_out = str(float_in / (np.pi))+"π"
return string_out
def format_secs_radians_label(float_in, cycle_length):
# https://stackoverflow.com/questions/21172228/python-matplotlib-polar-plots-with-angular-labels-in-radians
# Converts a float value in radians into a
# string representation of that float
cycle_ratio = float_in/(2*np.pi)
string_out = f'{np.round(cycle_ratio*cycle_length,2)}s\n({float_in / (np.pi)}π)'
return string_out
def getImage(path):
"""Utility function for reading image file"""
return OffsetImage(plt.imread(path), zoom=0.1)
def image_3d(ax, arr, label, xy, offset_x, offset_y, label_offset_y=0, zoom=0.05, pad=0):
""" Place an image (arr) as annotation at position xy
https://stackoverflow.com/questions/48180327/matplotlib-3d-scatter-plot-with-images-as-annotations
"""
im = OffsetImage(arr, zoom=zoom)
im.image.axes = ax
ab = AnnotationBbox(im, xy, xybox=(offset_x, offset_y),
xycoords='data', boxcoords="offset points",
pad=pad, arrowprops=dict(arrowstyle="->"))
ax.add_artist(ab)
offsetbox = TextArea(label, minimumdescent=False)
ab = AnnotationBbox(offsetbox, xy,
xybox=(offset_x, offset_y+label_offset_y),
xycoords='data',
boxcoords=("offset points"))
ax.add_artist(ab)
def plot_sorted_corr_mat(corr_mat, cluster_assignments, ax):
"""
Plot correlation matrix sorted by cluster assignments of nodes (must pass matplotlib axis)
"""
# Create sorting index from factor assignments
sort_indx = np.argsort(cluster_assignments)
sorted_vals = np.sort(cluster_assignments)
# Sort Distance Matrix
sortedmat = [[corr_mat[i][j] for j in sort_indx] for i in sort_indx]
# Plot Distance Matrix
c = ax.pcolormesh(sortedmat, cmap='coolwarm')
plt.colorbar(c, ax=ax)
# Plot rectangular patches along diagnol to indicate factor assignments
for i in np.unique(sorted_vals):
ind = np.where(sorted_vals == i)
mn = np.min(ind)
mx = np.max(ind)
sz=(mx-mn)+1
rect = Rectangle((mn,mn), sz, sz , linewidth=1,
edgecolor='black', facecolor='none')
ax.add_patch(rect)
def proj_3d(X, ax1, ax2):
""" From a 3D point in axes ax1,
calculate position in 2D in ax2
https://stackoverflow.com/questions/48180327/matplotlib-3d-scatter-plot-with-images-as-annotations
"""
x,y,z = X
x2, y2, _ = proj3d.proj_transform(x,y,z, ax1.get_proj())
return ax2.transData.inverted().transform(ax1.transData.transform((x2, y2)))
def shiftedColorMap(cmap, start=0, midpoint=0.5, stop=1.0, name='shiftedcmap'):
'''
Function to offset the "center" of a colormap. Useful for
data with a negative min and positive max and you want the
middle of the colormap's dynamic range to be at zero
Input
-----
cmap : The matplotlib colormap to be altered
start : Offset from lowest point in the colormap's range.
Defaults to 0.0 (no lower ofset). Should be between
0.0 and 1.0.
midpoint : The new center of the colormap. Defaults to
0.5 (no shift). Should be between 0.0 and 1.0. In
general, this should be 1 - vmax/(vmax + abs(vmin))
For example if your data range from -15.0 to +5.0 and
you want the center of the colormap at 0.0, `midpoint`
should be set to 1 - 5/(5 + 15)) or 0.75
stop : Offset from highets point in the colormap's range.
Defaults to 1.0 (no upper ofset). Should be between
0.0 and 1.0.
'''
cdict = {
'red': [],
'green': [],
'blue': [],
'alpha': []
}
# regular index to compute the colors
reg_index = np.linspace(start, stop, 257)
# shifted index to match the data
shift_index = np.hstack([
np.linspace(0, midpoint, 128, endpoint=False),
np.linspace(midpoint, 1.0, 129, endpoint=True)
])
for ri, si in zip(reg_index, shift_index):
r, g, b, a = cmap(ri)
cdict['red'].append((si, r, r))
cdict['green'].append((si, g, g))
cdict['blue'].append((si, b, b))
cdict['alpha'].append((si, a, a))
newcmap = LinearSegmentedColormap(name, cdict)
plt.register_cmap(cmap=newcmap)
return newcmap
def transition_matrix(transitions):
"""Generate markov transition from a 1-d sequence of state labels (numpy array)"""
#https://stackoverflow.com/questions/46657221/generating-markov-transition-matrix-in-python
n = 1 + np.int(np.nanmax(transitions)) #number of states
M = [[0]*n for _ in range(n)]
for (i,j) in zip(transitions,transitions[1:]):
if ~np.isnan(i) and ~np.isnan(j):
M[np.int(i)][np.int(j)] += 1
#now convert to probabilities:
for row in M:
s = sum(row)
if s > 0:
row[:] = [f/s for f in row]
return np.array(M)
def traveling_index(real_vec, imag_vec):
"""
Compute traveling index as the reciprocal of the condition number between the
real and imaginary component from cpca
"""
return 1/np.linalg.cond(np.vstack((real_vec, imag_vec)).T)
@np.vectorize
def twod_gauss(x, y, var=2):
"""Create two-dimensional gaussian with mean (mu) and variance parameters in the x- and y-plane"""
return np.exp(-((x - 0)**2 + (y - 0)**2)/(2*(var)**2))
def xcorr(x, y, maxlags=30):
"""Calculate cross-correlation from two time series (numpy array)"""
Nx = len(x)
if Nx != len(y):
raise ValueError('x and y must be equal length')
c = np.correlate(x, y, mode=2)
c /= np.sqrt(np.dot(x, x) * np.dot(y, y))
if maxlags is None:
maxlags = Nx - 1
if maxlags >= Nx or maxlags < 1:
raise ValueError('maglags must be None or strictly '
'positive < %d' % Nx)
lags = np.arange(-maxlags, maxlags + 1)
c = c[Nx - 1 - maxlags:Nx + maxlags]
max_r = c[np.argsort(np.abs(c))[-1]]
max_lag = lags[np.argsort(np.abs(c))[-1]]
return max_r, max_lag
# %% [markdown]
# ## Create ROY-BIG-BL brain colormap that matches HCP Workbench
# %%
colors_roy = ["cyan", "lime", "blueviolet", "mediumblue", "black", "red", "orange", "yellow"]
roy_big_bl = LinearSegmentedColormap.from_list("roy_big_bl", colors_roy)
# %% [markdown]
# # Figure 2 - Form and Properties of Three Dominant Spatiotemporal Patterns.
# %% [markdown]
# ## 1. Calculation of Complex Principal Component Duration and Time-Scale
# %%
pca_complex_res = pickle.load(open('results/pca_rest_complex_results.pkl', 'rb'))
cpca_recon = pickle.load(open('results/cpca_reconstruction_results.pkl', 'rb'))
# Unwrap temporal phase angles
comp0_phase = np.unwrap(np.angle(pca_complex_res['pca']['pc_scores'][:,0]))
comp1_phase = np.unwrap(np.angle(pca_complex_res['pca']['pc_scores'][:,1]))
comp2_phase = np.unwrap(np.angle(pca_complex_res['pca']['pc_scores'][:,2]))
# Calculate average duration of complete cycle from temporal phase angles
avg_cycle_comp0 = np.abs(comp0_phase[-1]-comp0_phase[0])/60000
avg_cycle_comp0 = (2*np.pi)/avg_cycle_comp0
avg_cycle_comp1 = np.abs(comp1_phase[-1]-comp1_phase[0])/60000
avg_cycle_comp1 = (2*np.pi)/avg_cycle_comp1
avg_cycle_comp2 = np.abs(comp2_phase[-1]-comp2_phase[0])/60000
avg_cycle_comp2 = (2*np.pi)/avg_cycle_comp2
# Calculate the range (in proportion of 2pi) of spatial phase angles
comp0_phase_weights = np.angle(pca_complex_res['pca']['Va'][0,:])
comp1_phase_weights = np.angle(pca_complex_res['pca']['Va'][1,:])
comp2_phase_weights = np.angle(pca_complex_res['pca']['Va'][2,:])
comp0_spatial_phase_range = (max(comp0_phase_weights) - min(comp0_phase_weights))/(2*np.pi)
comp1_spatial_phase_range = (max(comp1_phase_weights) - min(comp1_phase_weights))/(2*np.pi)
comp2_spatial_phase_range = (max(comp2_phase_weights) - min(comp2_phase_weights))/(2*np.pi)
# Calculate duration of lead-lag relationships in spatial phase angles
comp0_phase_weights_duration = comp0_spatial_phase_range * avg_cycle_comp0
comp1_phase_weights_duration = comp1_spatial_phase_range * avg_cycle_comp1
comp2_phase_weights_duration = comp2_spatial_phase_range * avg_cycle_comp2
# Derive time point units (in secs) for the spatial phase maps from the CPCA reconstruction (N = 30 bins)
comp0_phase_weights_units = (comp0_phase_weights_duration/30)*tr
comp1_phase_weights_units = (comp1_phase_weights_duration/30)*tr
comp2_phase_weights_units = (comp2_phase_weights_duration/30)*tr
# %% [markdown]
# ## 2. Figure
# %%
eigs = pickle.load(open('demo_files/pca_complex_eigenvalues.pkl', 'rb'))
exp_var = [eig/(n_vertices*2) for eig in eigs]
_, pca_comps, _ = pull_cifti_data(load_cifti('demo_files/pca_rest_complex_ang.dtseries.nii'))
zero_mask = np.std(pca_comps, axis=0) > 0
pca_comps = pca_comps[:, zero_mask]
hsv_cmap = plt.get_cmap('hsv')
fig = plt.figure(figsize=(18,20), constrained_layout=False)
crop_width_l = (0, 1000)
crop_width_r = (1300, 2100)
crop_height = (0,1053)
gspec = fig.add_gridspec(3,1, hspace=0.2, height_ratios=[0.33, 0.33, 0.33])
## Component 1
g_sub0 = gridspec.GridSpecFromSubplotSpec(2,2, subplot_spec=gspec[0], wspace=0,
width_ratios=[0.4, 0.6], height_ratios=[0.01,0.99])
title_ax = fig.add_subplot(g_sub0[0,0])
title_ax.set_title('A) First Complex Principal Component - Pattern One',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
# Compute travel index
real_vec = np.real(pca_complex_res['pca']['loadings'][0,:]).T
imag_vec = np.imag(pca_complex_res['pca']['loadings'][0,:]).T
t_index = traveling_index(real_vec, imag_vec)
title_ax = fig.add_subplot(g_sub0[0,1])
title_ax.set_title(f'Travel index:{np.round(t_index,2)}',
fontsize=14, loc='center')
title_ax.axis('off')
g_sub0_0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0[1:,0], hspace=0.1,
height_ratios = [0.6,0.4])
g_sub0_0_0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0_0[0], hspace=0.1,
height_ratios = [0.01,0.99])
g_sub0_0_1 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0_0[1], hspace=0.9,
height_ratios = [0.01,0.99])
title_ax = fig.add_subplot(g_sub0_0_0[0])
title_ax.set_title('Phase Delay Map',
fontsize=15, loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub0_0_1[0])
title_ax.set_title('Phase Delay Values',
fontsize=13, loc='left')
title_ax.axis('off')
ax = fig.add_subplot(g_sub0_0_0[1])
img = mpimg.imread('demo_files/pca_rest_complex_comp0_ang_hsv.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,(0,960)))
ax.axis('off')
box = ax.get_position()
box.y0 = box.y0 + 0.005
box.y1 = box.y1 + 0.005
ax.set_position(box)
xval = np.arange(-2.41, 2.6, 0.01)
yval = np.ones_like(xval)
hsv_cmap_shift0 = shiftedColorMap(hsv_cmap, midpoint=0.6)
norm0 = colors.Normalize(-2.42, 2.6)
ax = plt.subplot(g_sub0_0_1[1], polar=True)
ax.scatter(xval, yval, c=xval, s=300, cmap=hsv_cmap_shift0, norm=norm0, linewidths=0)
ax.set_yticks([])
convert_polar_xticks_to_radians_secs(ax, comp0_phase_weights_duration)
box = ax.get_position()
box.y0 = box.y0 + 0.005
box.y1 = box.y1 + 0.005
ax.set_position(box)
ax.tick_params(axis='both', which='major', pad=10)
g_sub0_1 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0[1:,1], hspace=0.05,
height_ratios = [0.01,0.99])
title_ax = fig.add_subplot(g_sub0_1[0])
title_ax.set_title('Reconstructed Time Points',
fontsize=15, loc='left')
title_ax.axis('off')
g_sub0_1_0 = gridspec.GridSpecFromSubplotSpec(2,3, subplot_spec=g_sub0_1[1], hspace=0.05,
wspace=0.05)
base_dir = 'demo_files/cpca_recon_pics'
t_samples = [0, 5, 9, 15, 20, 26]
for i, t in enumerate(t_samples):
ax = fig.add_subplot(g_sub0_1_0[i])
img = mpimg.imread(f'{base_dir}/cpca_recon_comp0_t{t}.png')
ax.set_title(f'{np.round(comp0_phase_weights_units * t,1)} s', fontsize=13)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
# Component 2
g_sub0 = gridspec.GridSpecFromSubplotSpec(2,2, subplot_spec=gspec[1], wspace=0,
width_ratios=[0.4, 0.6], height_ratios=[0.01,0.99])
title_ax = fig.add_subplot(g_sub0[0,0])
title_ax.set_title('B) Second Complex Principal Component - Pattern Two',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
# Compute travel index
real_vec = np.real(pca_complex_res['pca']['loadings'][1,:]).T
imag_vec = np.imag(pca_complex_res['pca']['loadings'][1,:]).T
t_index = traveling_index(real_vec, imag_vec)
title_ax = fig.add_subplot(g_sub0[0,1])
title_ax.set_title(f'Travel index:{np.round(t_index,2)}',
fontsize=14, loc='center')
title_ax.axis('off')
g_sub0_0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0[1:,0], hspace=0.1,
height_ratios = [0.6,0.4])
g_sub0_0_0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0_0[0], hspace=0.1,
height_ratios = [0.01,0.99])
g_sub0_0_1 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0_0[1], hspace=0.9,
height_ratios = [0.01,0.99])
title_ax = fig.add_subplot(g_sub0_0_0[0])
title_ax.set_title('Phase Delay Map',
fontsize=15, loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub0_0_1[0])
title_ax.set_title('Phase Delay Values',
fontsize=13, loc='left')
title_ax.axis('off')
ax = fig.add_subplot(g_sub0_0_0[1])
img = mpimg.imread('demo_files/pca_rest_complex_comp1_ang_hsv.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,(0,960)))
ax.axis('off')
box = ax.get_position()
box.y0 = box.y0 + 0.005
box.y1 = box.y1 + 0.005
ax.set_position(box)
xval = np.arange(-2.41, 2.6, 0.01)
yval = np.ones_like(xval)
hsv_cmap_shift0 = shiftedColorMap(hsv_cmap, midpoint=0.6)
norm0 = colors.Normalize(-2.42, 2.6)
ax = plt.subplot(g_sub0_0_1[1], polar=True)
ax.scatter(xval, yval, c=xval, s=300, cmap=hsv_cmap_shift0, norm=norm0, linewidths=0)
ax.set_yticks([])
convert_polar_xticks_to_radians_secs(ax, comp0_phase_weights_duration)
box = ax.get_position()
box.y0 = box.y0 + 0.005
box.y1 = box.y1 + 0.005
ax.set_position(box)
ax.tick_params(axis='both', which='major', pad=10)
g_sub0_1 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0[1:,1], hspace=0.05,
height_ratios = [0.01,0.99])
title_ax = fig.add_subplot(g_sub0_1[0])
title_ax.set_title('Reconstructed Time Points',
fontsize=15, loc='left')
title_ax.axis('off')
g_sub0_1_0 = gridspec.GridSpecFromSubplotSpec(2,3, subplot_spec=g_sub0_1[1], hspace=0.05,
wspace=0.05)
base_dir = 'demo_files/cpca_recon_pics'
t_samples = [0, 5, 10, 14, 19, 25]
for i, t in enumerate(t_samples):
ax = fig.add_subplot(g_sub0_1_0[i])
img = mpimg.imread(f'{base_dir}/cpca_recon_comp1_t{t}.png')
ax.set_title(f'{np.round(comp0_phase_weights_units * t,1)} s', fontsize=13)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
# Component 3
g_sub0 = gridspec.GridSpecFromSubplotSpec(2,2, subplot_spec=gspec[2], wspace=0,
width_ratios=[0.4, 0.6], height_ratios=[0.01,0.99])
title_ax = fig.add_subplot(g_sub0[0,0])
title_ax.set_title('C) Third Complex Principal Component - Pattern Three',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
# Compute travel index
real_vec = np.real(pca_complex_res['pca']['loadings'][2,:]).T
imag_vec = np.imag(pca_complex_res['pca']['loadings'][2,:]).T
t_index = traveling_index(real_vec, imag_vec)
title_ax = fig.add_subplot(g_sub0[0,1])
title_ax.set_title(f'Travel index:{np.round(t_index,2)}',
fontsize=14, loc='center')
title_ax.axis('off')
g_sub0_0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0[1:,0], hspace=0.1,
height_ratios = [0.6,0.4])
g_sub0_0_0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0_0[0], hspace=0.1,
height_ratios = [0.01,0.99])
g_sub0_0_1 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0_0[1], hspace=0.9,
height_ratios = [0.01,0.99])
title_ax = fig.add_subplot(g_sub0_0_0[0])
title_ax.set_title('Phase Delay Map',
fontsize=15, loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub0_0_1[0])
title_ax.set_title('Phase Delay Values',
fontsize=13, loc='left')
title_ax.axis('off')
ax = fig.add_subplot(g_sub0_0_0[1])
img = mpimg.imread('demo_files/pca_rest_complex_comp2_ang_hsv.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,(0,960)))
ax.axis('off')
box = ax.get_position()
box.y0 = box.y0 + 0.005
box.y1 = box.y1 + 0.005
ax.set_position(box)
xval = np.arange(-2.41, 2.6, 0.01)
yval = np.ones_like(xval)
hsv_cmap_shift0 = shiftedColorMap(hsv_cmap, midpoint=0.6)
norm0 = colors.Normalize(-2.42, 2.6)
ax = plt.subplot(g_sub0_0_1[1], polar=True)
ax.scatter(xval, yval, c=xval, s=300, cmap=hsv_cmap_shift0, norm=norm0, linewidths=0)
ax.set_yticks([])
convert_polar_xticks_to_radians_secs(ax, comp0_phase_weights_duration)
box = ax.get_position()
box.y0 = box.y0 + 0.005
box.y1 = box.y1 + 0.005
ax.set_position(box)
ax.tick_params(axis='both', which='major', pad=10)
g_sub0_1 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=g_sub0[1:,1], hspace=0.05,
height_ratios = [0.01,0.99])
title_ax = fig.add_subplot(g_sub0_1[0])
title_ax.set_title('Reconstructed Time Points',
fontsize=15, loc='left')
title_ax.axis('off')
g_sub0_1_0 = gridspec.GridSpecFromSubplotSpec(2,3, subplot_spec=g_sub0_1[1], hspace=0.05,
wspace=0.05)
base_dir = 'demo_files/cpca_recon_pics'
t_samples = [0, 5, 11, 15, 20, 25]
for i, t in enumerate(t_samples):
ax = fig.add_subplot(g_sub0_1_0[i])
img = mpimg.imread(f'{base_dir}/cpca_recon_comp2_t{t}.png')
ax.set_title(f'{np.round(comp0_phase_weights_units * t,1)} s', fontsize=13)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax.annotate('time',
xy=(0.28,0.75),
xytext=(0.24, 0.79),
xycoords='figure fraction',
textcoords='figure fraction',
fontweight='bold',
fontsize=12,
arrowprops=dict(facecolor='black', arrowstyle='<-', connectionstyle="angle3, angleA=0, angleB=80", alpha = 0.9, linewidth=3),
horizontalalignment='center',
verticalalignment='center',
bbox=dict(pad=5, facecolor="none", edgecolor="none")
)
plt.savefig('results/figures/cpca.eps', bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Figure Caption
# %% [markdown]
# Figure 2. Form and Properties of Three Dominant Spatiotemporal Patterns. Time-lag delay maps and reconstructed time points of the first three complex principal components. Time-lag delay maps represent the temporal ordering (in seconds) of cortical vertex BOLD time series within the spatiotemporal pattern. Time-lag delay maps describe a repeating or cyclical pattern expressed in radians (0 to 2) around a unit circle, where a phase value of 0 corresponds to the beginning of the spatiotemporal pattern, and 2 corresponds to the end of the spatiotemporal pattern. For clarity, radians are converted to temporal units (seconds) (see ‘Methods and Materials’). The values in the time-lag delay map correspond to the temporal delay (in seconds) between two cortical vertices, such that smaller values occur before larger values. Values are mapped to a cyclical color map to emphasize the cyclical temporal progression of each spatiotemporal pattern. To illustrate the temporal progression of the spatiotemporal patterns, six reconstructed time points are displayed for each pattern. A) The time-lag delay map (top) and reconstructed time points (bottom) of the first spatiotemporal pattern - ‘pattern one’. B) The time-lag delay map and reconstructed time points of the second spatiotemporal pattern - ‘pattern two’ C) The time-lag delay map and reconstructed time points of the third spatiotemporal pattern - ‘pattern three’.
# %% [markdown]
# # Figure 3 - Survey of Zero-lag FC Topographies.
# %%
cifti_fps = (
'demo_files/pca_rest.dtseries.nii',
'demo_files/eigenmap_p90.dtseries.nii',
'demo_files/fc_map_precuneus.dtseries.nii',
'demo_files/fc_map_sm.dtseries.nii',
'demo_files/fc_map_supramarginal.dtseries.nii',
'demo_files/pca_rest_varimax.dtseries.nii',
'demo_files/s_ica.dtseries.nii', 'demo_files/t_ica.dtseries.nii',
'demo_files/caps_precuneus_c2.dtseries.nii',
'demo_files/caps_sm_c2.dtseries.nii',
'demo_files/caps_supramarginal_c2.dtseries.nii',
'demo_files/hmm_mean_map.dtseries.nii'
)
labels_short = (
'Eigenmap 1', 'P Seed', 'SM Seed', 'SMG Seed',
'Varimax Comp 1', 'Varimax Comp 2', 'Varimax Comp 3', 'SICA Comp 1', 'SICA Comp 2',
'SICA Comp 3', 'TICA Comp 1', 'TICA Comp 2', 'TICA Comp 3', 'P CAP 1', 'P CAP 2',
'SM CAP 1', 'SM CAP 2', 'SMG CAP 1', 'SMG CAP 2', 'HMM State 1', 'HMM State 2',
'HMM State 3'
)
section_labels = (
'Laplacian Eigenmaps',
'Precuneus Seed Regression',
'Somatosensory Seed Regression',
'Supramarginal Seed Regression',
'Principal Component Analysis - Varimax Rotated',
'Spatial Independent Component Analysis',
'Temporal Independent Component Analysis',
'Precuneus Seed CAPS',
'SM Seed CAPs',
'Supramarginal Seed CAPS',
'GMM Hidden Markov Model'
)
## 1. Load All Maps
cifti_maps_all_orig = []
for fp in cifti_fps:
_, cifti_maps, n_time = pull_cifti_data(load_cifti(fp))
if any([label in fp for label in ['pca_rest', 'ica', 'hmm']]):
cifti_maps_all_orig.append(cifti_maps[:3, :])
elif 'eigenmap' in fp:
cifti_maps_all_orig.append(cifti_maps[0, :])
else:
cifti_maps_all_orig.append(cifti_maps[:2, :])
cifti_maps_all_orig = np.vstack(cifti_maps_all_orig)
zero_mask = np.std(cifti_maps_all_orig, axis=0) > 0
zero_mask_indx = np.where(zero_mask)[0]
cifti_maps_all = cifti_maps_all_orig[:, zero_mask].copy()
# Normalize
cifti_maps_all = zscore(cifti_maps_all.T)
## 2. Correlate all maps with first three principal components
component_corrs = np.corrcoef(cifti_maps_all.T)[3:,:3]
## 5. Calculate PCA explained variance
pca_ts = pickle.load(open('demo_files/pca_ts.pkl', 'rb'))
eigs = [np.var(pca_ts[:,i]) for i in range(pca_ts.shape[1])]
exp_var = [eig/n_vertices for eig in eigs]
## 6. Create Figure
crop_width_l = (0, 1000)
crop_width_r = (1300, 2100)
crop_height = (0,1053)
fig = plt.figure(figsize=(13,21), constrained_layout=False)
gspec = fig.add_gridspec(2,2, hspace=0.1, wspace=0.2,
width_ratios=[0.3, 0.7],
height_ratios=[0.8,0.2])
g_sub0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=gspec[:,0],
wspace=0, hspace=0.1,
height_ratios=[0.01,0.99])
g_sub1 = gridspec.GridSpecFromSubplotSpec(4,1, subplot_spec=gspec[0,1],
wspace=0, hspace=0.2,
height_ratios=[0.01,0.33,0.33,0.33])
g_sub2 = gridspec.GridSpecFromSubplotSpec(1,1, subplot_spec=gspec[1,1],
wspace=0)
title_ax = fig.add_subplot(g_sub0[0])
title_ax.set_title('A) Spatial Correlation with \n Principal Components',
fontsize=16, fontweight='bold', loc='left')
box = title_ax.get_position()
box.y0 = box.y0 - 0.02
box.y1 = box.y1 - 0.02
title_ax.set_position(box)
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub1[0])
title_ax.set_title('B) Principal Component Maps',
fontsize=16, fontweight='bold', loc='right')
box = title_ax.get_position()
box.y0 = box.y0 - 0.01
box.y1 = box.y1 - 0.01
title_ax.set_position(box)
title_ax.axis('off')
title_ax.axis('off')
ax = fig.add_subplot(g_sub2[0])
ax.set_aspect(0.009)
ax.plot(list(range(1,11)), eigs, '-o', color='black')
ax.set_xlabel('Component Number', fontsize=15)
ax.set_ylabel('Eigenvalue', fontsize=15)
x_adj = 0.4
y_adj = [-25,40,-30]
for i in range(3):
ax.text((i+1)+x_adj,eigs[i]+y_adj[i],
r'{}%'.format(np.round(exp_var[i]*100,1)),
fontsize=15, bbox=dict(facecolor='white', alpha=0.5))
ax.text(4.5, 940, 'Explained Variance', fontsize=15, bbox=dict(facecolor='white', alpha=0.5))
ax.set_title('C) Explained Variance Scree Plot',
fontsize=16, fontweight='bold', loc='left', pad=20)
box = ax.get_position()
box.x0 = box.x0 + 0.1
box.x1 = box.x1 + 0.1
ax.set_position(box)
pca_maps = ['demo_files/pca_rest_comp0.png',
'demo_files/pca_rest_comp1.png',
'demo_files/pca_rest_comp2.png']
weights_df = pd.DataFrame(component_corrs, columns=[f'Comp{i}' for i in range(3)], index=labels_short)
weights_df_abs = weights_df.abs()
df_list = []
indx_max = weights_df_abs.idxmax(axis=1)
for col in weights_df_abs.columns:
comp_weights = weights_df_abs.loc[indx_max==col]
comp_weights_sorted = comp_weights.sort_values(by=col, ascending=False)
df_list.append(comp_weights_sorted)
weights_df_abs_sorted = pd.concat(df_list)
ax = fig.add_subplot(g_sub0[1])
im = ax.imshow(weights_df_abs_sorted.values, aspect='auto')
ax.set_xticks(np.arange(3))
ax.set_yticks(np.arange(len(labels_short)))
ax.set_yticklabels(weights_df_abs_sorted.index, fontsize=15)
ax.set_xticklabels([f'PC{i+1}' for i in range(3)], fontsize=15,
fontweight='bold')
# Loop over data dimensions and create text annotations.
for i in range(weights_df_abs_sorted.shape[0]):
for j in range(weights_df_abs_sorted.shape[1]):
text = ax.text(j, i, round(weights_df_abs_sorted.iloc[i,j],2),
ha="center", va="center", color="ivory",
fontweight='bold', fontsize=13)
ax.tick_params(top=True, bottom=False, labeltop=True, labelbottom=False)
box = ax.get_position()
box.x0 = box.x0 + 0.1
box.x1 = box.x1 + 0.1
ax.set_position(box)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="10%", pad="20%")
cax.set_aspect(30)
cbar = fig.colorbar(im, cax=cax, orientation="vertical", aspect=10)
cbar.ax.tick_params(labelsize=14)
cbar.ax.set_title('Correlation', fontsize=15, loc='left')
ax1 = fig.add_subplot(g_sub1[1])
img = mpimg.imread(pca_maps[0])
ax1.set_title(f'Principal Component 1', fontsize=18)
ax1.imshow(cropImage(img, crop_width_l, crop_width_r, crop_height))
ax1.axis('off')
box = ax1.get_position()
box.x0 = box.x0 + 0.1
box.x1 = box.x1 + 0.1
ax1.set_position(box)
ax2 = fig.add_subplot(g_sub1[2])
img = mpimg.imread(pca_maps[1])
ax2.set_title(f'Principal Component 2', fontsize=18)
ax2.imshow(cropImage(img, crop_width_l, crop_width_r, crop_height))
ax2.axis('off')
box = ax2.get_position()
box.x0 = box.x0 + 0.1
box.x1 = box.x1 + 0.1
ax2.set_position(box)
ax3 = fig.add_subplot(g_sub1[3])
img = mpimg.imread(pca_maps[2])
ax3.set_title(f'Principal Component 3', fontsize=18)
ax3.imshow(cropImage(img, crop_width_l, crop_width_r, crop_height))
ax3.axis('off')
box = ax3.get_position()
box.x0 = box.x0 + 0.1
box.x1 = box.x1 + 0.1
ax3.set_position(box)
# plt.show()
plt.savefig('results/figures/FC_survey.eps', bbox_inches='tight')
# %% [markdown]
# ## Figure 3 Caption
# %% [markdown]
# Figure 3. Form and Properties of Three Fundamental Functional Connectivity Topographies. A) The spatial correlation (abs. value) between the first three principal component maps and each FC topography displayed as a table. The color of each cell in the table is shaded from light yellow (strong correlation) to dark blue (weak correlation). All FC topographies in our survey exhibited strong spatial correlations (Pearson’s correlation) with one (or two) of the first three principal components. B) The first three principal component spatial maps. C) The scree plot that displays the explained variance in cortical time series for each successive principal component. The scree plot indicates a clear elbow after the third principal component, indicating a ‘diminishing return’ in explained variance of extracting more components. (P=precuneus, SM= somatosensory; SMG=supramarginal gyrus; Clus=cluster; Comp=component; PC = Principal Component).
# %% [markdown]
# ## Permutation Spin Tests of FC Topographies
# %% [markdown]
# Compute permuation spin tests for each FC topograhy and most correlated principal component
# %%
l_sphere = load_surf_mesh('templates/sphere_left.gii')
r_sphere = load_surf_mesh('templates/sphere_right.gii')
# Get top principal component corr
sim_pairs = weights_df_abs_sorted.idxmax(axis=1)
pc_sim_pair = list(zip(sim_pairs,sim_pairs.index))
# Number of permutations per test
n_rand = 1000
# Run permutation spin test per map
labels = list(labels_short)
labels = ['Comp0', 'Comp1', 'Comp2'] + labels
perm_r_all = []
orig_r = []
for pair in pc_sim_pair:
print(pair)
# Index cifti maps
pc_indx = labels.index(pair[0])
map_indx = labels.index(pair[1])
pc_cifti = cifti_maps_all_orig[pc_indx, :]
map_cifti = cifti_maps_all_orig[map_indx, :]
orig_r.append(np.corrcoef(pc_cifti, map_cifti)[0,1])
# Split cifti map into left and right hemispheres
map_cifti_L, map_cifti_R = map_cifti[:n_vert_L], map_cifti[n_vert_L:]
# Initialize permutation test
sp = SpinPermutations(n_rep=n_rand, random_state=0)
sp.fit(l_sphere.coordinates, r_sphere.coordinates)
# randomize
map_rotated = np.hstack(sp.randomize(map_cifti_L, map_cifti_R))
# Calculate permutation distribution
perm_r = []
for i in range(n_rand):
perm_r.append(np.corrcoef(pc_cifti, map_rotated[i,:])[0,1])
perm_r_all.append(perm_r)
# Calculate p-value per map
pv_all = []
for pair, perm_r, r_obs in zip(pc_sim_pair, perm_r_all, orig_r):
# Include observed value as part of permutation distribution
pv = ( (np.abs(perm_r) >= np.abs(r_obs)).sum() + 1 )/(n_rand + 1)
pv_all.append((pair, pv))
# %% [markdown]
# # Movie 1 - Visualization of Spatiotemporal Patterns
# %%
%%HTML
<video controls autoplay loop>
<source
src="demo_files/time_lag_structures.mp4"
type="video/mp4"
</video>
# %% [markdown]
# ## Movie 1 Caption
# %% [markdown]
# Visualization of Spatiotemporal Patterns. Temporal reconstruction of all three spatiotemporal patterns displayed as movies in the following order - pattern one, pattern two, and pattern three. The time points are equally-spaced samples (N=30) of the spatiotemporal patterns. The seconds since the beginning of the spatiotemporal pattern are displayed in the top left. In the bottom of the panel, the time points of the spatiotemporal pattern are displayed in three-dimensional principal component space (Figure 2). Two-dimensional slices of the three principal component space (see Figure 3) are displayed as the three 2-dimensional plots. The progression of time points in the principal component space is illustrated by a cyclical color map (light to dark to light). The movement of the spatiotemporal pattern through this space is illustrated by a moving red dot from time point-to-time point in synchronization with the temporal reconstruction in the movie.
# %% [markdown]
# # Movie 2. Dynamic Visualization of the Quasiperiodic Pattern, Pattern One, and Global Signal.
# %%
%%HTML
<video controls autoplay loop>
<source
src="demo_files/qpp_comparison.mp4"
type="video/mp4"
</video>
# %% [markdown]
# ## Movie 2 Caption
# %% [markdown]
# Movie 2. Dynamic Visualization of the Quasiperiodic Pattern, Pattern One and Global Signal. The 30 time points (TR=0.72s) of the QPP, pattern one, and peak-average global signal displayed as a movie (in that order). The time index of each sequence is displayed in the top left. The time points of pattern one are equally-spaced phase samples (N=30) of the time point reconstruction (see above). The time points of the QPP are derived from the spatiotemporal template computed from the repeated-template-averaging procedure. The global signal visualization concatenates the left and right windows (w=15TRs) of the global signal peak-average. The time points of the global signal visualization begin at TR=-15, corresponding to 15 TRs pre-peak, and proceed to TR=15, corresponding to 15TRs post-peak.
# %% [markdown]
# # **Figure 4 - Similar Propagation Patterns between Average Latency Structure and Pattern One**
# %%
_, pca_comps, _ = pull_cifti_data(load_cifti('demo_files/pca_rest_complex_ang.dtseries.nii'))
_, lag_proj, _ = pull_cifti_data(load_cifti('demo_files/lag_projection.dtseries.nii'))
_, phase_circ_mean, _ = pull_cifti_data(load_cifti('demo_files/phase_circular_mean.dtseries.nii'))
zero_mask = np.std(pca_comps, axis=0) > 0
pca_comps = pca_comps[:, zero_mask]
lag_proj = lag_proj[:, zero_mask][0,:]
phase_circ_mean = phase_circ_mean[:, zero_mask]
corr_circ_lag = np.corrcoef(phase_circ_mean, lag_proj)[0,1]
corr_circ_p1 = np.corrcoef(phase_circ_mean, pca_comps[0,:])[0,1]
fig = plt.figure(figsize=(20,10), constrained_layout=False)
crop_width_l = (0, 1000)
crop_width_r = (1300, 2100)
crop_height = (0,1053)
gspec = fig.add_gridspec(1,3, wspace=0.2, width_ratios=[0.33,0.33,0.33])
ax = fig.add_subplot(gspec[0])
img = mpimg.imread('demo_files/lag_projection.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax.set_title('Lag Projection Map', fontsize=16, fontweight='bold')
ax = fig.add_subplot(gspec[1])
img = mpimg.imread('demo_files/phase_circular_mean.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax.set_title('Circular Average of Complex Correlations', fontsize=16, fontweight='bold')
ax = fig.add_subplot(gspec[2])
img = mpimg.imread('demo_files/pca_rest_complex_comp0_ang.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax.set_title('Pattern One Phase Delay Map', fontsize=16, fontweight='bold')
fig.text(0.355, 0.52, f"r = {round(corr_circ_lag, 2)}", va="center", fontsize=20)
fig.text(0.635, 0.52, f"r = {round(corr_circ_p1, 2)}", va="center", fontsize=20)
fig.text(0.205, 0.38, r"TR (0.72s)", va="center", fontsize=16)
fig.text(0.48, 0.38, r"$\Theta$ (radian)", va="center", fontsize=16)
fig.text(0.755, 0.38, r"$\Theta$ (radian)", va="center", fontsize=16)
plt.savefig('results/figures/cpca_lag_proj.eps', bbox_inches='tight')
plt.show()
# %%
np.corrcoef(phase_circ_mean, lag_proj)[0,1]
# %%
np.corrcoef(phase_circ_mean, pca_comps[0,:])[0,1]
# %% [markdown]
# ## Figure 4 Caption
# %% [markdown]
# Figure 4. Similar Time-lag Dynamics between Pattern One and Lag Projection. Comparison between the phase delay map of pattern one (left) and the average lag projection (right). As in Figure 2, the pattern one phase delay map represents the phase delay (in radians) of cortical BOLD time series. The lag projection map represents the average time-lag delay (in seconds) between each vertex of the cortex. The spatial correlation between the pattern one phase delay map and lag projection is r = 0.81, indicating a strong similarity in time-lag dynamics.
# %% [markdown]
# # Figure 5 - The Task-Positive/Task-Negative Pattern, Primary Gradient, and Pattern Two Describe the Same Spatiotemporal Pattern.
# %%
gs_signal = pickle.load(open('demo_files/gs_results.pkl', 'rb'))
comp_ts = pickle.load(open('demo_files/pca_complex_ts.pkl', 'rb'))
comp0_ts = np.real(comp_ts[:,0])
gs_comp0_corr = np.corrcoef(gs_signal, comp0_ts)[0,1]
eigenmap_indx = np.arange(0,100,10)
eigenmaps = []
for indx in eigenmap_indx:
_, eigenmap, _ = pull_cifti_data(load_cifti(f'demo_files/eigenmap_p{indx}.dtseries.nii'))
eigenmaps.append(eigenmap[0,:])
eigenmaps_array = np.array(eigenmaps)
zero_mask = np.std(eigenmaps_array, axis=0) > 0
eigenmaps_array = eigenmaps_array[:,zero_mask].copy()
_, pca_maps, _ = pull_cifti_data(load_cifti(f'demo_files/pca_rest.dtseries.nii'))
pca_maps = pca_maps[:2, zero_mask].copy()
pca_eigenmap_corr = np.corrcoef(pca_maps, eigenmaps_array)[:2, 2:]
fig = plt.figure(figsize=(22,18), constrained_layout=False)
crop_width_l = (0, 1000)
crop_width_r = (1300, 2100)
crop_height = (0,1053)
gspec = fig.add_gridspec(6,1, hspace=0.4, wspace=0,
height_ratios=[0.02,0.33,0.02,0.33,0.02,0.33])
title_ax = fig.add_subplot(gspec[0,:])
title_ax.set_title('A) Pattern Two, TP/TN Pattern & Primary Functional Connectivity Gradient',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
g_sub0 = gridspec.GridSpecFromSubplotSpec(1,3, subplot_spec=gspec[1],
wspace=0)
ax = fig.add_subplot(g_sub0[0])
img = mpimg.imread('demo_files/pca_rest_comp1_flip.png')
ax.set_title(f'Second Principal Component \n (sign flipped)', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax = fig.add_subplot(g_sub0[1])
img = mpimg.imread('demo_files/fc_map_precuneus_gs_nonsymmetric.png')
ax.set_title(f'Task-Positive/Task-Negative Pattern \n (precuneus seed)', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax = fig.add_subplot(g_sub0[2])
img = mpimg.imread('demo_files/eigenmap_p0_comp0.png')
ax.set_title(f'Primary Functional Connectivity Gradient \n (no threshold)', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
title_ax = fig.add_subplot(gspec[2,:])
title_ax.set_title('B) PCA & cPCA of Time-Point Centered Data',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
g_sub1 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=gspec[3], wspace=0, width_ratios=[0.62,0.38])
g_sub1_0 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=g_sub1[0], wspace=0)
g_sub1_1 = gridspec.GridSpecFromSubplotSpec(1,1, subplot_spec=g_sub1[1])
ax = fig.add_subplot(g_sub1_0[0])
img = mpimg.imread('demo_files/pca_rest_comp0.png')
ax.set_title(f'First Principal Component \n (original)', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax = fig.add_subplot(g_sub1_0[1])
img = mpimg.imread('demo_files/pca_rest_comp0_tmode_flip.png')
ax.set_title(f'First Principal Component \n (time-point centered)', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
hsv_cmap = plt.get_cmap('hsv')
ax = fig.add_subplot(g_sub1_1[0])
img = mpimg.imread('demo_files/pca_rest_complex_tmode_comp0_ang_hsv.png')
ax.set_title('First Complex Principal Component - \n Time-Lag Delay Map \n (time-point centered)', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
title_ax = fig.add_subplot(gspec[4,:])
title_ax.set_title('C) Threshold Effect on Primary Functional Connectivity Gradient',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
g_sub2 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=gspec[5], wspace=0.1, width_ratios=[0.25,0.75])
g_sub2_0 = gridspec.GridSpecFromSubplotSpec(1,1, subplot_spec=g_sub2[0], wspace=0)
g_sub2_1 = gridspec.GridSpecFromSubplotSpec(1,5, subplot_spec=g_sub2[1], wspace=0)
ax = fig.add_subplot(g_sub2_0[0])
ax.plot(eigenmap_indx, np.abs(pca_eigenmap_corr[0,:]), marker='o',color='b', label='PC1')
ax.plot(eigenmap_indx, np.abs(pca_eigenmap_corr[1,:]), marker='o',color='r', label='PC2')
ax.legend()
ax.set_xticks(eigenmap_indx)
ax.set_title('Spatial Correlation of Laplacian Eigenmap with \n First Two Principal Components', fontsize=15, pad=15)
ax.set_xlabel('Percentile Threshold of Affinity Matrix', fontsize=15)
ax.set_ylabel('Spatial Correlation (abs)', fontsize=15)
crop_width2 = (35, 1200)
crop_height2 = (5,20)
sub_title_ax = fig.add_subplot(g_sub2_1[0,:])
sub_title_ax.set_title('Laplacian Eigenmap by Percentile Threshold (Left Hemisphere)',
fontsize=16, loc='left')
sub_title_ax.axis('off')
box = sub_title_ax.get_position()
box.y0 = box.y0 + 0.015
box.y1 = box.y1 + 0.015
sub_title_ax.set_position(box)
eigenmap_indx2 = np.arange(0,100,20)
for i, indx in enumerate(eigenmap_indx2):
ax = fig.add_subplot(g_sub2_1[i])
img = mpimg.imread(f'demo_files/eigenmap_p{indx}_comp0.png')
ax.set_title(f'{indx}%', fontsize=14)
ax.imshow(cropImage_single(img,crop_width2,crop_height2))
ax.axis('off')
# plt.show()
plt.savefig('results/figures/fpn_to_dmn.eps', bbox_inches='tight')
# %% [markdown]
# ## Figure 5 Caption
# %% [markdown]
# Figure 5. The Task-Positive/Task-Negative Pattern, Primary Gradient, and Pattern Two Describe the Same Spatiotemporal Pattern. A) From left to right, pattern two, task-positive/task-negative (TP/TN) pattern, and the PG represented by the spatial weights of the second principal component from PCA (sign flipped for consistency), seed-based correlation map (precuneus seed), and first Laplacian eigenmap with no thresholding of the affinity matrix, respectively. As can be observed visually, similar spatial patterns are produced from all three analyses - pattern two:TP/TN (r =0.96) and pattern two:PG (r = 0.83). B) From left to right, the first principal component from non-time-centered BOLD time courses (i.e. pattern one), the first principal component of time-centered BOLD time courses, and the first complex principal component time-lag delay map from time-centered BOLD time courses. As can be observed visually, time-point centering of BOLD time courses replaces the original unipolar first principal component (left; pattern one) with a bipolar (anti-correlated) principal component (middle) that resembles pattern two. In the same manner, the first complex principal component of CPCA of time-centered BOLD time courses (right) exhibits a time-lag map resembling the time-lag map of pattern two (Figure 2). C) The effect of functional connectivity (FC) matrix percentile thresholding on the resulting spatial weights of the PG, computed as the first eigenmap of the Laplacian Eigenmap algorithm (only the left hemisphere presented for space). At zero to low-thresholding of the FC matrix, the first Laplacian Eigenmap resembles pattern two (PC2). As the threshold is raised, the spatial weights of vertices within the FPN, DMN and SMLV become more uniform, and the spatial weights of the vertices within the FPN fall to zero. At higher thresholds this results in an Eigenmap that resembles pattern one.
# %% [markdown]
# # Figure 6 - Comparison of Original and Reconstructed Functional Connectivity Matrices.
# %%
pca_res = pickle.load(open('results/pca_rest_complex_results.pkl', 'rb'))['pca']
orig_corr = pickle.load(open('results/fc_matrix_results.pkl', 'rb'))
U = pca_res['U'][:,:3]
s = np.diag(pca_res['s'][:3])
Va = pca_res['Va'][:3,:]
recon_ts = np.real(U @ s @ Va)
recon_corr = np.corrcoef(recon_ts.T)
orig_tril = orig_corr[np.tril_indices(orig_corr.shape[0], k=1)]
recon_tril = recon_corr[np.tril_indices(recon_corr.shape[0],k=1)]
orig_recon_corr = np.corrcoef(orig_tril, recon_tril)[0,1]
orig_corr_adj = orig_corr.copy()
recon_corr_adj = recon_corr.copy()
np.fill_diagonal(orig_corr_adj, 0)
np.fill_diagonal(recon_corr_adj, 0)
# # Louvain Community detection - original and reconstructed corr matrices
labels_orig, _ = community_louvain(orig_corr_adj, B='negative_asym', seed=0)
labels_recon, _ = community_louvain(recon_corr_adj, B='negative_asym', seed=0)
# Calculate normalized mutual information
_, norm_mni = partition_distance(labels_orig, labels_recon)
# # Write module assignments to cifti file for display
# # This cannot be done without access to data, so it is commented out
# # ex_subj_file = ['data/rest/proc_3_surf_resamp/200513_LR1_rest_smooth_filt_smooth_filt_resamp.R.func.gii',
# # 'data/rest/proc_3_surf_resamp/200513_LR1_rest_smooth_filt_smooth_filt_resamp.L.func.gii']
# # hdr = load_gifti(ex_subj_file)
# # write_to_gifti(labels_orig[np.newaxis, :]+1, hdr, 'louvain_modules', zero_mask)
fig = plt.figure(figsize=(25,5), constrained_layout=False)
# Define 1 by 3 overall grid
gspec = fig.add_gridspec(1,2, wspace=0.05, width_ratios=[0.65,0.35])
# Plot original correlation matrix with community partition
gspec0 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=gspec[0], wspace=0.5)
ax0 = fig.add_subplot(gspec0[0])
plot_sorted_corr_mat(orig_corr, labels_orig, ax0)
ax0.set_title('Original FC Matrix', fontsize=16, pad=8, fontweight='bold')
ax0.set_xlabel('vertex', fontsize=12)
ax0.set_ylabel('vertex', fontsize=12)
ax0.text(500, 2100, 'Module 1', fontweight='bold', fontsize=12)
ax0.text(2450, 3700, 'Module 2', fontweight='bold', fontsize=12)
ax0.text(3850, 4900, 'Module 3', fontweight='bold', fontsize=12)
# Plot reconstructed correlation matrix with community partition
ax1 = fig.add_subplot(gspec0[1])
plot_sorted_corr_mat(recon_corr, labels_recon, ax1)
ax1.set_title('Reconstructed FC Matrix', fontsize=16, pad=8, fontweight='bold')
ax1.set_xlabel('vertex', fontsize=12)
ax1.set_ylabel('vertex', fontsize=12)
fig.text(0.33, 0.8, f"r = {round(orig_recon_corr, 2)}", va="center", fontsize=20)
fig.text(0.325, 0.65, f"NMI = {round(norm_mni, 2)}", va="center", fontsize=20)
# Display module assignments on cortical surface
g_sub = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=gspec[1], wspace=0,
width_ratios=[0.7, 0.3])
ax2 = fig.add_subplot(g_sub[0])
img = mpimg.imread('demo_files/louvain_modules.png')
ax2.set_title(f'Module Partition', fontsize=16, fontweight='bold')
ax2.imshow(img)
ax2.axis('off')
# Display module labels
g_sub_sub = gridspec.GridSpecFromSubplotSpec(3,1, subplot_spec=g_sub[1], hspace = 0.1)
color_module = ['black', 'red', 'yellow']
for i, (gsub, clr) in enumerate(zip(g_sub_sub, color_module)):
sub_ax = fig.add_subplot(gsub)
sub_ax.text(0.4,0.5, f'Module {i+1}', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.2,0.4),0.05,0.3,linewidth=1,
edgecolor='none',facecolor=clr,
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
plt.savefig('results/figures/corr_compare.png', dpi=200)
# plt.show()
# %% [markdown]
# ### Figure 6 Caption
# %% [markdown]
# Figure 6. The Network Structure of Functional Connectivity Is Explained by the Three Fundamental Spatiotemporal Patterns. Comparison of the correlation matrix of cortical BOLD time courses (left) with the correlation matrix of reconstructed cortical BOLD time courses (right) derived from the three spatiotemporal patterns, and the module assignments of each vertex (bottom). The rows and columns of the original and reconstructed correlation matrix are sorted and outlined (in black) according to the modular structure estimated from the Louvain modularity algorithm. The algorithm identified three primary modules in the SMLV, DMN and FPN. Despite a higher mean value of correlations in the reconstructed correlation matrix, the pattern of correlations between the two correlation matrices is highly similar (r = 0.77). Further, the modular structure of the original correlation matrix exhibits a high degree of similarity with the modular structure of the reconstructed correlation matrix (MNI = 0.73).
# %% [markdown]
# # Supplementary Figures
# %% [markdown]
# # Supplementary Results A - Spatiotemporal Patterns Consist of Steady States and Propagation Events That Repeat Across Patterns.
# %%
V = pickle.load(open('demo_files/pca_eigen.pkl', 'rb')) # X = USV
V = zscore(V.T).T
_, cpca_comp0, n_time = pull_cifti_data(load_cifti('results/cpca_comp0_recon.dtseries.nii'))
_, cpca_comp1, n_time = pull_cifti_data(load_cifti('results/cpca_comp1_recon.dtseries.nii'))
_, cpca_comp2, n_time = pull_cifti_data(load_cifti('results/cpca_comp2_recon.dtseries.nii'))
zero_mask = np.std(cpca_comp0, axis=0) > 0
cpca_comp0 = zscore(cpca_comp0[:, zero_mask].copy().T).T
cpca_comp1 = zscore(cpca_comp1[:, zero_mask].copy().T).T
cpca_comp2 = zscore(cpca_comp2[:, zero_mask].copy().T).T
cpca_comp0_proj = cpca_comp0 @ V.T
cpca_comp1_proj = cpca_comp1 @ V.T
cpca_comp2_proj = cpca_comp2 @ V.T
fig = plt.figure(figsize=(20,25), constrained_layout=False)
crop_width_l = (0, 1000)
crop_width_r = (1300, 2100)
crop_height = (0,1053)
gspec = fig.add_gridspec(4,1, hspace=0.3, wspace=0,
height_ratios=[0.25,0.3,0.3,0.4])
g_sub0 = gridspec.GridSpecFromSubplotSpec(2,3, subplot_spec=gspec[0], wspace=0,
hspace=0.2, height_ratios=[0.01,0.99])
title_ax = fig.add_subplot(g_sub0[0,:])
title_ax.set_title('A) Principal Components',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
ax = fig.add_subplot(g_sub0[1,0])
img = mpimg.imread('demo_files/pca_rest_comp0.png')
ax.set_title('PC 1', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax = fig.add_subplot(g_sub0[1,1])
img = mpimg.imread('demo_files/pca_rest_comp1.png')
ax.set_title('PC 2', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax = fig.add_subplot(g_sub0[1,2])
img = mpimg.imread('demo_files/pca_rest_comp2.png')
ax.set_title('PC 3', fontsize=16)
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
g_sub1 = gridspec.GridSpecFromSubplotSpec(2,3, subplot_spec=gspec[1], wspace=0.2,
hspace=0, height_ratios=[0.01,0.99],
width_ratios=[0.3,0.3,0.4])
title_ax = fig.add_subplot(g_sub1[0,:])
title_ax.set_title('B) Spatiotemporal Patterns in Principal Component Space',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
t = np.arange(30)
ax = fig.add_subplot(g_sub1[1,0])
ax.set_xlabel('PC 1', fontsize=16, fontweight='bold')
ax.set_ylabel('PC 2', fontsize=16, fontweight='bold', labelpad=-20)
ax.scatter(cpca_comp0_proj[:,0], cpca_comp0_proj[:,1], c=t,
cmap='Blues', s=100, alpha=0.8)
ax.plot(cpca_comp0_proj[:,0], cpca_comp0_proj[:,1], alpha=0.5, color='b',
label='Pattern One')
ax.scatter(cpca_comp1_proj[:,0], cpca_comp1_proj[:,1], c=t,
cmap='Greens', s=100, alpha=0.5)
ax.plot(cpca_comp1_proj[:,0], cpca_comp1_proj[:,1], alpha=0.5, color='g',
label='Pattern Two')
ax.scatter(cpca_comp2_proj[:,0], cpca_comp2_proj[:,1], c=t,
cmap='Reds', s=100, alpha=0.5)
ax.plot(cpca_comp2_proj[:,0], cpca_comp2_proj[:,1], alpha=0.5, color='r',
label='Pattern Three')
ax.legend()
ax = fig.add_subplot(g_sub1[1,1])
ax.set_xlabel('PC 1', fontsize=16, fontweight='bold')
ax.set_ylabel('PC 3', fontsize=16, fontweight='bold', labelpad=-20)
ax.scatter(cpca_comp0_proj[:,0], cpca_comp0_proj[:,2], c=t,
cmap='Blues', s=100, alpha=0.8)
ax.plot(cpca_comp0_proj[:,0], cpca_comp0_proj[:,2], alpha=0.5, color='b')
ax.scatter(cpca_comp1_proj[:,0], cpca_comp1_proj[:,2], c=t,
cmap='Greens', s=100, alpha=0.8)
ax.plot(cpca_comp1_proj[:,0], cpca_comp1_proj[:,2], alpha=0.5, color='g')
ax.scatter(cpca_comp2_proj[:,0], cpca_comp2_proj[:,2], c=t,
cmap='Reds', s=100, alpha=0.8)
ax.plot(cpca_comp2_proj[:,0], cpca_comp2_proj[:,2], alpha=0.5, color='r')
ax = fig.add_subplot(g_sub1[1,2])
ax.set_xlabel('PC 2', fontsize=16, fontweight='bold')
ax.set_ylabel('PC 3', fontsize=16, fontweight='bold', labelpad=-20)
scatter1 = ax.scatter(cpca_comp0_proj[:,1], cpca_comp0_proj[:,2], c=t,
cmap='Blues', s=100, alpha=0.8)
ax.plot(cpca_comp0_proj[:,1], cpca_comp0_proj[:,2], alpha=0.5, color='b')
scatter2 = ax.scatter(cpca_comp1_proj[:,1], cpca_comp1_proj[:,2], c=t,
cmap='Greens', s=100, alpha=0.8)
ax.plot(cpca_comp1_proj[:,1], cpca_comp1_proj[:,2], alpha=0.5, color='g')
scatter3 = ax.scatter(cpca_comp2_proj[:,1], cpca_comp2_proj[:,2], c=t,
cmap='Reds', s=100, alpha=0.8)
ax.plot(cpca_comp2_proj[:,1], cpca_comp2_proj[:,2], alpha=0.5, color='r')
cbar = plt.colorbar(scatter3, shrink=0.7, pad=0, ax=ax, fraction=0.1)
cbar.ax.set_title('Pattern One TR', loc='left', rotation=50, fontsize=13)
cbar.ax.tick_params(labelsize=13)
cbar = plt.colorbar(scatter2, shrink=0.7, pad=0, ax=ax, fraction=0.1)
cbar.ax.set_title('Pattern Two TR', loc='left', rotation=50, fontsize=13)
cbar.set_ticks([])
cbar = plt.colorbar(scatter1, shrink=0.7, pad=0.01, ax=ax, fraction=0.1)
cbar.ax.set_title('Pattern Three TR', loc='left', rotation=50, fontsize=13)
cbar.set_ticks([])
cpca_comp_all = np.vstack([zscore(cpca_comp0.T).T, zscore(cpca_comp1.T).T, zscore(cpca_comp2.T).T])
cpca_comp_all_proj = np.vstack([cpca_comp0_proj, cpca_comp1_proj, cpca_comp2_proj])
kmeans = KMeans(n_clusters=6, n_init=10, random_state=0)
kmeans.fit(cpca_comp_all)
g_sub2 = gridspec.GridSpecFromSubplotSpec(2,3, subplot_spec=gspec[2], wspace=0.2,
hspace=0.1, width_ratios=[0.3,0.3,0.4],
height_ratios=[0.01,0.99])
title_ax = fig.add_subplot(g_sub2[0,:])
title_ax.set_title('C) Recurring Pattern Clusters in Spatiotemporal Patterns',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
ax0 = fig.add_subplot(g_sub2[1,0])
ax1 = fig.add_subplot(g_sub2[1,1])
ax2 = fig.add_subplot(g_sub2[1,2])
ax0.set_xlabel('PC 1', fontsize=16, fontweight='bold')
ax0.set_ylabel('PC 2', fontsize=16, fontweight='bold', labelpad=-20)
ax1.set_xlabel('PC 1', fontsize=16, fontweight='bold')
ax1.set_ylabel('PC 3', fontsize=16, fontweight='bold', labelpad=-20)
ax2.set_xlabel('PC 2', fontsize=16, fontweight='bold')
ax2.set_ylabel('PC 3', fontsize=16, fontweight='bold', labelpad=-20)
c_colors = [plt.cm.tab20(i) for i in [14,15,16,17,18,19]]
for c_indx, label in enumerate(np.unique(kmeans.labels_)):
indx = np.where(kmeans.labels_==label)
ax0.scatter(cpca_comp_all_proj[indx,0], cpca_comp_all_proj[indx,1],
color=c_colors[c_indx], s=100, label=f'Cluster {c_indx+1}')
ax1.scatter(cpca_comp_all_proj[indx,0], cpca_comp_all_proj[indx,2],
color=c_colors[c_indx], s=100)
ax2.scatter(cpca_comp_all_proj[indx,1], cpca_comp_all_proj[indx,2],
color=c_colors[c_indx], s=100)
ax0.plot(cpca_comp0_proj[:,0], cpca_comp0_proj[:,1], alpha=0.2, color='b')
ax0.plot(cpca_comp1_proj[:,0], cpca_comp1_proj[:,1], alpha=0.2, color='g')
ax0.plot(cpca_comp2_proj[:,0], cpca_comp2_proj[:,1], alpha=0.2, color='r')
ax0.legend()
ax1.plot(cpca_comp0_proj[:,0], cpca_comp0_proj[:,2], alpha=0.2, color='b')
ax1.plot(cpca_comp1_proj[:,0], cpca_comp1_proj[:,2], alpha=0.2, color='g')
ax1.plot(cpca_comp2_proj[:,0], cpca_comp2_proj[:,2], alpha=0.2, color='r')
ax2.plot(cpca_comp0_proj[:,1], cpca_comp0_proj[:,2], alpha=0.2, color='b')
ax2.plot(cpca_comp1_proj[:,1], cpca_comp1_proj[:,2], alpha=0.2, color='g')
ax2.plot(cpca_comp2_proj[:,1], cpca_comp2_proj[:,2], alpha=0.2, color='r')
labels_c0 = kmeans.labels_[:30]
labels_c1 = kmeans.labels_[30:60]
labels_c2 = kmeans.labels_[60:]
labels_df = pd.DataFrame([labels_c0, labels_c1, labels_c2]).T
cmap = colors.ListedColormap(c_colors)
bounds=[0,1,2,3,4,5]
norm = colors.BoundaryNorm(bounds, cmap.N)
ax2_divider = make_axes_locatable(ax2)
sub_ax2 = ax2_divider.append_axes("right", size="30%", pad="10%")
sub_ax2.imshow(labels_df, cmap=cmap)
sub_ax2.set_aspect(0.5)
sub_ax2.set_ylabel('TR', fontsize=14, labelpad=-2)
sub_ax2.set_title('D) Clusters by TR', fontsize=18, fontweight='bold', pad=25)
sub_ax2.set_xticks([0,1,2])
sub_ax2.set_xticklabels(['Pattern One', 'Pattern Two', 'Pattern Three'], fontsize=12,
rotation=50, ha='center')
## Write cluster centroids to cifti files
# This cannot be done without access to data, so it is commented out
# ex_subj_file = ['data/rest/proc_3_surf_resamp/200513_LR1_rest_smooth_filt_smooth_filt_resamp.R.func.gii',
# 'data/rest/proc_3_surf_resamp/200513_LR1_rest_smooth_filt_smooth_filt_resamp.L.func.gii']
# hdr = load_gifti(ex_subj_file)
# write_to_gifti(kmeans.cluster_centers_, hdr, 'cpca_kmeans_N6', zero_mask)
g_sub3 = gridspec.GridSpecFromSubplotSpec(3,3, subplot_spec=gspec[3], wspace=0.1,
hspace=0, height_ratios=[0.01,0.49,0.49])
title_ax = fig.add_subplot(g_sub3[0,:])
title_ax.set_title('E) Recurring Pattern Cluster Maps',
fontsize=18, fontweight='bold', loc='left')
title_ax.axis('off')
ax = fig.add_subplot(g_sub3[1,0])
img = mpimg.imread('demo_files/cpca_recon_cluster0.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax_divider = make_axes_locatable(ax)
sub_ax = ax_divider.append_axes("top", size="10%", pad="2%")
sub_ax.text(0.4,0, 'Cluster 1', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.7,0),0.05,1,linewidth=1,
edgecolor='none',facecolor=c_colors[0],
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
ax = fig.add_subplot(g_sub3[1,1])
img = mpimg.imread('demo_files/cpca_recon_cluster1.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax_divider = make_axes_locatable(ax)
sub_ax = ax_divider.append_axes("top", size="10%", pad="2%")
sub_ax.text(0.4,0, 'Cluster 2', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.7,0),0.05,1,linewidth=1,
edgecolor='none',facecolor=c_colors[1],
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
ax = fig.add_subplot(g_sub3[1,2])
img = mpimg.imread('demo_files/cpca_recon_cluster2.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax_divider = make_axes_locatable(ax)
sub_ax = ax_divider.append_axes("top", size="10%", pad="2%")
sub_ax.text(0.4,0, 'Cluster 3', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.7,0),0.05,1,linewidth=1,
edgecolor='none',facecolor=c_colors[2],
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
ax = fig.add_subplot(g_sub3[2,0])
img = mpimg.imread('demo_files/cpca_recon_cluster3.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax_divider = make_axes_locatable(ax)
sub_ax = ax_divider.append_axes("top", size="10%", pad="2%")
sub_ax.text(0.4,0, 'Cluster 4', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.7,0),0.05,1,linewidth=1,
edgecolor='none',facecolor=c_colors[3],
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
ax = fig.add_subplot(g_sub3[2,1])
img = mpimg.imread('demo_files/cpca_recon_cluster4.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax_divider = make_axes_locatable(ax)
sub_ax = ax_divider.append_axes("top", size="10%", pad="2%")
sub_ax.text(0.4,0, 'Cluster 5', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.7,0),0.05,1,linewidth=1,
edgecolor='none',facecolor=c_colors[4],
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
ax = fig.add_subplot(g_sub3[2,2])
img = mpimg.imread('demo_files/cpca_recon_cluster5.png')
ax.imshow(cropImage(img,crop_width_l,crop_width_r,crop_height))
ax.axis('off')
ax_divider = make_axes_locatable(ax)
sub_ax = ax_divider.append_axes("top", size="10%", pad="2%")
sub_ax.text(0.4,0, 'Cluster 6', fontsize=16)
# add a fancy box
fancybox = FancyBboxPatch((0.7,-0.1),0.05,1,linewidth=1,
edgecolor='none',facecolor=c_colors[5],
boxstyle=BoxStyle("Round", pad=0.02, rounding_size=0.015))
sub_ax.add_patch(fancybox)
sub_ax.axis('off')
# plt.show()
plt.savefig('results/figures/cpca_dynamics.eps', bbox_inches='tight')
# %% [markdown]
# ## Figure Caption
# %% [markdown]
# Supplementary Figure A. Spatiotemporal Patterns Consist of Steady States and Propagation Events That Repeat Across Patterns. (PC = Principal Component) Illustration of the progression of BOLD activity over time in each spatiotemporal pattern. A) The spatial weights for the first three principal components from PCA (Figure 1B). For visualization of temporal dynamics, the reconstructed time points (N=30) from each spatiotemporal pattern were projected onto the 3-dimensional embedding space formed by the first three principal components. B) Two-dimensional slices of each spatiotemporal pattern in the 3-dimensional principal component space - PC1-PC2, PC1-PC3, and PC2-PC3 spaces. The time points of patterns one, two and three are displayed as blue, green and red points, respectively. Consecutive time points of each spatiotemporal pattern are linked by lines. The time points of each spatiotemporal pattern are colored from light to darker to visualize the progression of time (N=30). The score of each time point on a given principal component is proportional to the Pearson correlation coefficient between the BOLD activity at that time point with the spatial weights of the principal component. Examination of the movement of time points within the 3-dimensional space provides information regarding the temporal dynamics of the spatiotemporal pattern. C) The same two-dimensional slices of each spatiotemporal pattern in the 3-dimensional principal component space colored according to their cluster assignment by a k-means clustering algorithm. K-means clustering was used to identify recurring spatial patterns of BOLD activity across time points of the three spatiotemporal patterns. Six clusters were estimated. D) The cluster assignments (color) by time (y-axis) of each spatiotemporal pattern (x-axis). Note, that the same cluster assignment can occur across more than one spatiotemporal pattern. E) The cluster centroids from the k-means clustering algorithm, corresponding to the average spatial pattern of BOLD activity for the time points that belong to that cluster. Note, the cluster centroids of the first two clusters are mean-centered versions of the original unimodal (all-positive or all-negative) steady-state of pattern one, as z-score normalization of the time-points across vertices was performed beforehand.
# %% [markdown]
# # Supplementary Figure A - Low-Dimensional Latent FC Topograhies
# %%
cifti_fps = (
'demo_files/pca_rest.dtseries.nii',
'demo_files/eigenmap_p90.dtseries.nii',
'demo_files/pca_rest_varimax.dtseries.nii',
'demo_files/s_ica.dtseries.nii', 'demo_files/t_ica.dtseries.nii'
)
fps = (
'demo_files/pca_rest_comp0.png', 'demo_files/pca_rest_comp1.png', 'demo_files/pca_rest_comp2.png',
'demo_files/eigenmap_p90_comp0.png',
'demo_files/pca_rest_varimax_comp0.png', 'demo_files/pca_rest_varimax_comp1.png',
'demo_files/pca_rest_varimax_comp2.png', 'demo_files/spatial_ica_comp0.png', 'demo_files/spatial_ica_comp1.png',
'demo_files/spatial_ica_comp2.png', 'demo_files/temporal_ica_comp0.png', 'demo_files/temporal_ica_comp1.png',
'demo_files/temporal_ica_comp2.png'
)
labels_short = (
'PCA Comp 1', 'PCA Comp 2', 'PCA Comp 3',
'Eigenmap 1', 'Varimax Comp 1', 'Varimax Comp 2', 'Varimax Comp 3',
'SICA Comp 1', 'SICA Comp 2', 'SICA Comp 3', 'TICA Comp 1', 'TICA Comp 2',
'TICA Comp 3'
)
pca_ts = pickle.load(open('demo_files/pca_ts.pkl', 'rb'))[:,:3]
varimax_ts = pickle.load(open('demo_files/varimax_ts.pkl', 'rb'))
tica_ts = pickle.load(open('demo_files/tica_ts.pkl', 'rb'))
sica_ts = pickle.load(open('demo_files/sica_ts.pkl', 'rb'))
all_ts = [pca_ts, varimax_ts, sica_ts, tica_ts]
corr_time = np.corrcoef(np.hstack(all_ts).T)
## 1. Load All Maps
cifti_maps_all = []
for fp in cifti_fps:
_, cifti_maps, n_time = pull_cifti_data(load_cifti(fp))
if any([label in fp for label in ['pca_rest', 'ica']]):
cifti_maps_all.append(cifti_maps[:3, :])
elif 'eigenmap' in fp:
cifti_maps_all.append(cifti_maps[0, :])
else:
cifti_maps_all.append(cifti_maps[:2, :])
cifti_maps_all = np.vstack(cifti_maps_all)
zero_mask = np.std(cifti_maps_all, axis=0) > 0
zero_mask_indx = np.where(zero_mask)[0]
cifti_maps_all = cifti_maps_all[:, zero_mask].copy()
# # Normalize
corr_maps = np.corrcoef(cifti_maps_all)
# cifti_maps_all = zscore(cifti_maps_all.T)
fig = plt.figure(figsize=(16,22), constrained_layout=False)
# gspec = fig.add_gridspec(5,3, hspace=0.2, wspace=0, height_ratios=[0.85,0.01,0.15])
gspec = fig.add_gridspec(9,3, hspace=0.2, wspace=0, height_ratios=[0.01,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15])
title_ax = fig.add_subplot(gspec[0,:])
title_ax.set_title('A) Low-Dimensional FC Topograhies',
fontsize=16, fontweight='bold', loc='left')
title_ax.axis('off')
# title_ax = fig.add_subplot(gspec[1,:])
# title_ax.set_title('B) Spatial and Temporal Correlations with Principal Components',
# fontsize=16, fontweight='bold', loc='left')
# title_ax.axis('off')
# g_sub0 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=gspec[0], hspace=0, height_ratios=[0.73,0.27])
# g_sub0_0 = gridspec.GridSpecFromSubplotSpec(2,4, subplot_spec=g_sub0[0], hspace=0, wspace=0)
# g_sub0_1 = gridspec.GridSpecFromSubplotSpec(1,3, subplot_spec=g_sub0[1], wspace=0)
crop_width = (35, 150)
crop_height = (5,15)
# indx=3
# for i in range(3):
# for j in range(4):
# if indx < 11:
# ax = fig.add_subplot(g_sub0_0[i,j])
# elif indx < 14:
# ax = fig.add_subplot(g_sub0_1[j])
# else:
# break
# img = mpimg.imread(fps[indx])
# ax.set_title(labels_short[indx], fontsize=16)
# ax.imshow(cropImage(img,crop_width,crop_height))
# ax.axis('off')
# indx+=1
indx=0
grid_indices = ([1,0], [1,1], [1,2], [2,0], [3,0], [4,0], [5,0], [6,0],
[6,1], [6,2], [7,0], [7,1], [7,2])
for g_indx in grid_indices:
ax = fig.add_subplot(gspec[g_indx[0], g_indx[1]])
img = mpimg.imread(fps[indx])
ax.set_title(labels_short[indx], fontsize=16)
ax.imshow(cropImage_single(img,crop_width,crop_height))
ax.axis('off')
indx+=1
# g_sub1 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=gspec[2])
ax = fig.add_subplot(gspec[3:6,1])
im = ax.imshow(np.abs(corr_maps[:3,3:].T), aspect=0.3, cmap='coolwarm', vmin=0, vmax=0.9)
ax.set_xticks(np.arange(3))
ax.set_yticks(np.arange(len(labels_short[3:])))
ax.set_yticklabels(labels_short[3:], fontweight='bold')
ax.set_xticklabels(labels_short[:3], rotation=30, fontweight='bold')
ax.xaxis.tick_top()
ax.set_title('B) Spatial Correlation with PCs', fontsize=16, fontweight='bold', loc='center')
ax.set_aspect(1)
box = ax.get_position()
box.x0 = box.x0 - 0.01
box.x1 = box.x1 - 0.01
box.y0 = box.y0 + 0.04
box.y1 = box.y1 + 0.04
ax.set_position(box)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="10%", pad=0.1)
cbar = plt.colorbar(im, orientation='vertical', cax=cax)
ax = fig.add_subplot(gspec[3:6,2])
im = ax.imshow(np.abs(corr_time[:3,3:].T), aspect=0.25, cmap='coolwarm', vmin=0, vmax=0.9)
ax.set_xticks(np.arange(3))
ax.set_yticks(np.arange(len(labels_short[4:])))
ax.set_yticklabels(labels_short[4:], fontweight='bold')
ax.set_xticklabels(labels_short[:3], rotation=30, fontweight='bold')
ax.xaxis.tick_top()
ax.set_title('C) Temporal Correlations with PCs', fontsize=16, fontweight='bold', loc='center')
ax.set_aspect(1)
box = ax.get_position()
box.y0 = box.y0 + 0.04
box.y1 = box.y1 + 0.04
ax.set_position(box)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="10%", pad=0.1)
cbar = plt.colorbar(im, orientation='vertical', cax=cax)
plt.savefig('results/figures/supplement_latentFC.eps')
plt.show()
# %% [markdown]
# ### Supplementary Figure A Caption
# %% [markdown]
# (SICA=Spatial ICA; TICA = Temporal ICA). The spatial weights of components from PCA (N=3), Laplacian Eigenmaps (N=1), varimax rotation of principal components (N=3), spatial ICA (N=3) and temporal ICA (N=3). The temporal and spatial correlations (absolute value) between the components of dimension-reduction analyses and the first three principal components are shown in the middle of the plot. Note, due to the nature of the Laplacian Eigenmap algorithm as a non-linear manifold learning algorithm, time courses cannot be extracted for their components. As illustrated in the spatial and temporal correlations table, the dimension-reduction analyses are largely consistent in their spatial topographies and temporal dynamics with the first three principal components.
# %% [markdown]
# # Supplementary Figure B - Seed-Based Topographies
# %%
fps = (
['demo_files/fc_map_sm.png', 'demo_files/fc_map_sm_gs.png'],
['demo_files/fc_map_precuneus.png', 'demo_files/fc_map_precuneus_gs.png'],
['demo_files/fc_map_supramarginal.png', 'demo_files/fc_map_supramarginal_gs.png'],
['demo_files/caps_sm_cluster0_c2.png', 'demo_files/caps_sm_cluster1_c2.png'],
['demo_files/caps_precuneus_cluster0_c2.png', 'demo_files/caps_precuneus_cluster1_c2.png'],
['demo_files/caps_supramarginal_cluster0_c2.png', 'demo_files/caps_supramarginal_cluster1_c2.png'],
['demo_files/caps_sm_norm_cluster0_c2.png', 'demo_files/caps_sm_norm_cluster1_c2.png'],
['demo_files/caps_precuneus_norm_cluster0_c2.png', 'demo_files/caps_precuneus_norm_cluster1_c2.png'],
['demo_files/caps_supramarginal_norm_cluster0_c2.png', 'demo_files/caps_supramarginal_norm_cluster1_c2.png']
)
pca_complex_res = pickle.load(open('results/pca_rest_complex_results.pkl', 'rb'))['pca']
comp_ts = np.real(pca_complex_res['pc_scores'][:,:3])
crop_width_L = (35, 1200)
crop_width_R = (1100, 150)
crop_height = (5,20)
fig = plt.figure(figsize=(20,28), constrained_layout=False)
gspec = fig.add_gridspec(2,1, hspace=0.05, wspace=0,
height_ratios=[0.6,0.4])
g_sub0 = gridspec.GridSpecFromSubplotSpec(9,6, hspace=0.2, wspace=0, subplot_spec=gspec[0],
height_ratios=[0.001,0.02,0.3,
0.001,0.02,0.3,
0.001,0.02,0.3])
g_sub1 = gridspec.GridSpecFromSubplotSpec(4,6, hspace=0.3, wspace=0, subplot_spec=gspec[1],
height_ratios=[0.01,0.4,
0.1,0.4])
title_ax = fig.add_subplot(g_sub0[0,:])
title_ax.set_title('A) Seed-Based Regression Maps',
fontsize=16, fontweight='bold', loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub0[3,:])
title_ax.set_title('B) Co-activation Pattern Clusters (N=2)',
fontsize=16, fontweight='bold', loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub0[6,:])
title_ax.set_title('C) Time-point Normalized - Co-actvation Pattern Clusters (N=2)',
fontsize=16, fontweight='bold', loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub1[0,:])
title_ax.set_title('D) Overlap in Suprathreshold Time Points',
fontsize=16, fontweight='bold', loc='left')
title_ax.axis('off')
title_ax = fig.add_subplot(g_sub1[2,:])
title_ax.text(0, 0,'E) Correlation between Suprathreshold Time Points and Three Time-lag Structures',
fontsize=16, fontweight='bold')
title_ax.axis('off')
axis_rows = [1,4,7]
for row in axis_rows:
sub_title = fig.add_subplot(g_sub0[row,:2])
sub_title.text(0.2, 0, 'Somatosensory Cortex Seed',
fontsize=14, fontweight='bold')
sub_title.axis('off')
sub_title = fig.add_subplot(g_sub0[row,2:4])
sub_title.text(0.3, 0, 'Precuneus Seed',
fontsize=14, fontweight='bold')
sub_title.axis('off')
sub_title = fig.add_subplot(g_sub0[row,4:6])
sub_title.text(0.2, 0, 'Supramarginal Gryus Seed',
fontsize=14, fontweight='bold')
sub_title.axis('off')
indx = 0
for fp_list in fps[:3]:
for fp in fp_list:
if indx % 2 == 0:
crop_w = crop_width_L
label='Original'
else:
crop_w = crop_width_R
label='Global Signal Regressed'
ax = fig.add_subplot(g_sub0[2,indx])
img = mpimg.imread(fp)
ax.set_title(label, fontsize=14)
ax.imshow(cropImage_single(img,crop_w,crop_height))
ax.axis('off')
indx+=1
indx = 0
for fp_list in fps[3:6]:
for fp in fp_list:
if indx % 2 == 0:
crop_w = crop_width_L
label='Cluster 1'
else:
crop_w = crop_width_R
label='Cluster 2'
ax = fig.add_subplot(g_sub0[5,indx])
img = mpimg.imread(fp)
ax.set_title(label, fontsize=14)
ax.imshow(cropImage_single(img,crop_w,crop_height))
ax.axis('off')
indx+=1
indx = 0
for fp_list in fps[6:]:
for fp in fp_list:
if indx % 2 == 0:
crop_w = crop_width_L
label='Cluster 1'
else:
crop_w = crop_width_R
label='Cluster 2'
ax = fig.add_subplot(g_sub0[8,indx])
img = mpimg.imread(fp)
ax.set_title(label, fontsize=14)
ax.imshow(cropImage_single(img,crop_w,crop_height))
ax.axis('off')
indx+=1
# Load and create CAP time series
seeds = ['sm', 'precuneus', 'supramarginal']
seed_labels = ['SM', 'P', 'SMG']
clus_ts = {}
clus_ts_norm = {}
for seed, seed_label in zip(seeds, seed_labels):
clus_ts[seed_label] = {}
clus_ts_norm[seed_label] = {}
caps_res = pickle.load(open(f'results/caps_{seed}_c2_results.pkl', 'rb'))
caps_res_norm = pickle.load(open(f'results/caps_{seed}_norm_c2_results.pkl', 'rb'))
ts_indx = caps_res[2]
clus_indx = caps_res[1]
ts_indx_norm = caps_res_norm[2]
clus_indx_norm = caps_res_norm[1]
for clus in [0,1]:
clus_ts_indx = ts_indx[clus_indx==clus]
clus_ts_indx_norm = ts_indx_norm[clus_indx_norm==clus]
ts_tmp = np.zeros(n_ts)
ts_tmp_norm = np.zeros(n_ts)
ts_tmp[clus_ts_indx] = 1
ts_tmp_norm[clus_ts_indx_norm] = 1
clus_ts[seed_label][f'C{clus+1}'] = ts_tmp
clus_ts_norm[seed_label][f'C{clus+1}'] = ts_tmp_norm
# Create dataframe of CAP time series
all_ts = []
all_ts_norm = []
all_ts_labels = []
for seed in seed_labels:
for clus in ['C1', 'C2']:
all_ts.append(clus_ts[seed][clus])
all_ts_norm.append(clus_ts_norm[seed][clus])
all_ts_labels.append(seed + '_' + clus)
all_ts = pd.DataFrame(np.array(all_ts).T, columns=all_ts_labels)
all_ts_norm = pd.DataFrame(np.array(all_ts_norm).T, columns=all_ts_labels)
# Calculate overlap in CAP time series w/ jaccard similarity
sim_mat = np.zeros((6,6))
sim_mat_norm = np.zeros((6,6))
sim_mat_lag = np.zeros((6,6))
sim_mat_lag_norm = np.zeros((6,6))
for x in range(6):
for y in range(x,6):
x_ts = all_ts.iloc[:, [x]]; x_ts_norm = all_ts_norm.iloc[:, [x]]
y_ts = all_ts.iloc[:, [y]]; y_ts_norm = all_ts_norm.iloc[:, [y]]
lags = list(range(-30,31))
lag_jaccard = [1 - cdist(x_ts.shift(i).values.T, y_ts.values.T, 'jaccard')[0]
for i in lags]
lag_jaccard_norm = [1 - cdist(x_ts_norm.shift(i).values.T, y_ts_norm.values.T, 'jaccard')[0]
for i in lags]
max_jaccard = np.max(lag_jaccard); max_jaccard_norm = np.max(lag_jaccard_norm)
max_jaccard_lag = lags[np.argmax(lag_jaccard)]
max_jaccard_lag_norm = lags[np.argmax(lag_jaccard_norm)]
sim_mat[x,y] = max_jaccard; sim_mat_norm[x,y] = max_jaccard_norm
sim_mat_lag[x,y] = max_jaccard_lag; sim_mat_lag_norm[x,y] = max_jaccard_lag_norm
i_lower = np.tril_indices(6, -1)
sim_mat[i_lower] = sim_mat.T[i_lower]
sim_mat_norm[i_lower] = sim_mat_norm.T[i_lower]
sim_mat_lag[i_lower] = sim_mat_lag.T[i_lower]
sim_mat_lag_norm[i_lower] = sim_mat_lag_norm.T[i_lower]
# Sort jaccard similarity matrix
labels = [1, 0, 1, 0, 0, 1]
sort_indx = np.argsort(labels)
sorted_vals = np.sort(labels)
labels_sorted = [all_ts_labels[i] for i in sort_indx]
sortedmat = [[sim_mat[i,j] for j in sort_indx] for i in sort_indx]
sim_mat_sorted = pd.DataFrame(sortedmat, columns = labels_sorted, index=labels_sorted)
sortedmat = [[sim_mat_norm[i,j] for j in sort_indx] for i in sort_indx]
sim_mat_sorted_norm = pd.DataFrame(sortedmat, columns = labels_sorted, index=labels_sorted)
g_sub1_0 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=g_sub1[1, :], wspace=0.5)
ax1 = fig.add_subplot(g_sub1_0[0])
ax2 = fig.add_subplot(g_sub1_0[1])
im1 = ax1.imshow(sim_mat_sorted, vmin=0, vmax=0.25)
ax1.set_xticks(np.arange(6))
ax1.set_yticks(np.arange(6))
ax1.set_yticklabels(labels_sorted, fontsize=13, fontweight='bold')
ax1.set_xticklabels(labels_sorted, fontsize=13, rotation=55, fontweight='bold')
ax1.set_title('D1) Cluster Jaccard Similarity', fontsize=14,
fontweight='bold', loc='center', pad=7)
plt.colorbar(im1, ax=ax1)
im2 = ax2.imshow(sim_mat_sorted_norm, vmin=0, vmax=0.25)
ax2.set_xticks(np.arange(6))
ax2.set_yticks(np.arange(6))
ax2.set_yticklabels(labels_sorted, fontsize=13, fontweight='bold')
ax2.set_xticklabels(labels_sorted, fontsize=13, rotation=55, fontweight='bold')
ax2.set_title('D2) Cluster Jaccard Similarity - Normalized', fontsize=14,
fontweight='bold', loc='center', pad=7)
plt.colorbar(im2, ax=ax2)
g_sub1_1 = gridspec.GridSpecFromSubplotSpec(1,2, subplot_spec=g_sub1[3, :], wspace=0.5)
ax1 = fig.add_subplot(g_sub1_1[0])
ax2 = fig.add_subplot(g_sub1_1[1])
# Calculate correlation between CAP time series and time-lag structures
corr_mat = np.zeros((3,6))
corr_mat_norm = np.zeros((3,6))
corr_mat_lag = np.zeros((3,6))
corr_mat_lag_norm = np.zeros((3,6))
for x in range(3):
for y in range(6):
pc_ts = pd.Series(comp_ts[:, x])
cap_ts = all_ts.iloc[:, y]; cap_ts_norm = all_ts_norm.iloc[:, y]
lags = list(range(-30,31))
lag_corr = np.abs([pc_ts.corr(cap_ts.shift(i)) for i in lags])
lag_corr_norm = np.abs([pc_ts.corr(cap_ts_norm.shift(i)) for i in lags])
max_corr = np.max(lag_corr); max_corr_norm = np.max(lag_corr_norm)
max_corr_lag = lags[np.argmax(lag_corr)]
max_corr_lag_norm = lags[np.argmax(lag_corr_norm)]
corr_mat[x,y] = max_corr; corr_mat_norm[x,y] = max_corr_norm
corr_mat_lag[x,y] = max_corr_lag; corr_mat_lag_norm[x,y] = max_corr_lag_norm
pc_labels = ['SMLV-to-FPN', 'FPN-to-DMN', 'FPN-to-SMLV']
im1 = ax1.imshow(corr_mat, vmin=0, vmax=0.4)
cax = plt.colorbar(im1, ax=ax1)
cax.ax.set_title('Correlation \n (abs. value)')
ax1.set_xticks(np.arange(6))
ax1.set_yticks(np.arange(3))
ax1.set_yticklabels(pc_labels, fontsize=13, fontweight='bold')
ax1.set_xticklabels(all_ts_labels, fontsize=13, rotation=55, fontweight='bold')
ax1.set_title('E1) Correlations b/w CAPs and \n Time-lag Structures', fontsize=14,
fontweight='bold', loc='center', pad=7)
ax1.set_aspect(0.5)
im2 = ax2.imshow(corr_mat_norm, vmin=0, vmax=0.4)
cax = plt.colorbar(im2, ax=ax2)
cax.ax.set_title('Correlation \n (abs. value)')
ax2.set_xticks(np.arange(6))
ax2.set_yticks(np.arange(3))
ax2.set_yticklabels(pc_labels, fontsize=13, fontweight='bold')
ax2.set_xticklabels(all_ts_labels, fontsize=13, rotation=55, fontweight='bold')
ax2.set_title('E2) Correlations b/w Normalized CAPs and \n Time-lag Structures', fontsize=14,
fontweight='bold', loc='center', pad=7)
ax2.set_aspect(0.5)
plt.savefig('results/figures/supplement_seeds.eps')
plt.show()
# %% [markdown]
# ### Supplementary Figure B Caption
# %% [markdown]
# (SM = somatosensory cortex; P=Precuneus; SMG=Supramarginal Gyrus). Spatial topographies of seed-based regression maps and CAP centroids from somatosensory (SM), precuneus and supramarginal gyrus seeds. A) Seed-based regression maps with (left hemisphere) and without global signal regression (right hemisphere) for SM, precuneus and supramarginal gyrus seeds. B) CAP cluster centroids (N=2) from k-means clustering of non-normalized (i.e. not z-scored) suprathreshold time points from SM, precuneus and supramarginal seeds. C) CAP cluster centroids (N=2) of the same suprathreshold time points with normalization (i.e. z-scored) before input to the k-means clustering algorithm. D1) Temporal overlap between binary time courses (see main text) of the two CAPs from each seed using the Jaccard similarity (Jaccard index). The Jaccard similarity between two CAP binary time courses varies from 0 to 1, and reflects the ratio of overlapping onset time points (=1) to the total number of time points (N=60,000). D2) Temporal overlap between CAP binary time courses from the normalized solutions of each seed analysis. E) Temporal correlation between the beginning phase time course of the three time-lag structures (SMLV-to-FPN, FPN-to-DMN and FPN-to-SMLV) and the CAP binary time courses for the non-normalized (E1) and normalized (E) solutions.
# %% [markdown]
# # Supplementary Figure D - Scree Plot from Complex Principal Component Analysis.
# %%
fig, ax = plt.subplots(figsize=(7,7))
eigs = pickle.load(open('demo_files/pca_complex_eigenvalues.pkl', 'rb'))
exp_var = [eig/(n_vertices*2) for eig in eigs]
ax.plot(list(range(1,11)), eigs, '-o', color='black')
ax.set_xlabel('Component Number', fontsize=13)
ax.set_ylabel('Eigenvalue', fontsize=13)
x_adj = 0.3
y_adj = [-20,50,-30]
for i in range(3):
ax.text((i+1)+x_adj,eigs[i]+y_adj[i],
r'{}%'.format(np.round(exp_var[i]*100,1)),
fontsize=13, bbox=dict(facecolor='white', alpha=0.5))
ax.text(6, 2000, 'Explained Variance', fontsize=13, bbox=dict(facecolor='white', alpha=0.5))
ax.set_title('Complex PCA Scree Plot', fontweight='bold', fontsize=17)
plt.savefig('results/figures/supplementaryC_screeplot.eps', bbox_inches='tight')
plt.show()
# %% [markdown]
# The eigenvalue by component number plot (i.e. scree plot) used to determine the number of components to extract. There are clear elbows in the plot after one and three components, indicating a preferred solution of one or three principal components (three were chosen).
# %% [markdown]
# # Supplementary Figure F. Principal Component and Functional Connectivity Gradient Topographies
# %%
fps = [['demo_files/pca_rest_comp0.png', 'demo_files/pca_rest_comp1.png', 'demo_files/pca_rest_comp2.png'],
['demo_files/pca_rest_gs_comp0.png', 'demo_files/pca_rest_gs_comp1.png', 'demo_files/pca_rest_gs_comp2.png'],
['demo_files/pca_rest_comp0_tmode.png', 'demo_files/pca_rest_comp1_tmode.png', 'demo_files/pca_rest_comp2_tmode.png'],
['demo_files/eigenmap_p0_comp0.png', 'demo_files/eigenmap_p0_comp1.png', 'demo_files/diffusion_emb_comp2.png']]
labels = [['Component 1', 'Component 2', 'Component 3'],
['Component 1', 'Component 2', 'Component 3'],
['Component 1', 'Component 2', 'Component 3'],
['Eigenmap 1', 'Eigenmap 2', 'Eigenmap 3']]
section_labels = ['Principal Component Analysis',
'Principal Component Analysis - Global Signal Removed',
'Principal Component Analysis - Time-Point Centered',
'Laplacian Eigenmaps - Manifold Learning']
fig = plt.figure(figsize=(20,20), constrained_layout=False)
gspec = fig.add_gridspec(8,3, hspace=0.05, wspace=0,
height_ratios=[0.01,0.2,0.01,0.2,0.01,0.2,0.01,0.2])
title_inds = [0,2,4,6]
for title_indx, label in zip(title_inds, section_labels):
title_ax = fig.add_subplot(gspec[title_indx,:])
title_ax.set_title(label,fontsize=16, fontweight='bold', loc='left')
title_ax.axis('off')
img_inds = [1,3,5,7]
for fp_sec, label_sec, img_indx in zip(fps, labels, img_inds):
for i in range(3):
ax = fig.add_subplot(gspec[img_indx,i])
img = mpimg.imread(fp_sec[i])
ax.set_title(label_sec[i], fontsize=15)
ax.imshow(img)
ax.axis('off')
fig.set_facecolor('w')
plt.savefig('results/figures/supplementaryD_pcagradients.eps')
plt.show()
# %% [markdown]
# ### Supplementary Figure D Caption
# %% [markdown]
# Displayed are the FC topography spatial weights from PCA, PCA on global-signal regressed data, PCA on time-point centered data, and Laplacian Eigenmaps. Note, we observed that the eigenmaps were highly positively skewed. To make the negative values of the eigenmaps more visible the colormap is made non-symmetric. The first and second eigenmaps match the second and third principal component from PCA. The first principal component is missing from the LE, global-signal regressed, and time-point centered PCA solutions.
# %% [markdown]
# # Supplementary Figure E. Comparison of Lag Projections With and Without Global Signal Regression.
# %%
fps = ['demo_files/lag_projection.png', 'demo_files/lag_projection_gs.png']
labels = ['Lag Projection - Without Global Signal Regression',
'Lag Projection - With Global Signal Regression']
fig, axs = plt.subplots(figsize=(15, 15) , nrows=1, ncols=2)
img = mpimg.imread(fps[0])
axs[0].imshow(img)
axs[0].set_title(labels[0])
axs[0].axis('off')
img = mpimg.imread(fps[1])
axs[1].imshow(img)
axs[1].set_title(labels[1])
axs[1].axis('off')
fig.set_facecolor('w')
plt.savefig('results/figures/supplementaryG_lag_projection.eps', bbox_inches='tight')
# plt.show()
# %% [markdown]
# ### Figure E Caption
# %% [markdown]
# Lag projections with and without global signal regression as a preprocessing step. Values on each cortical map represent the average time-delay between each cortical vertex and all others. Time-delay values are colored from light green/blue (earlier in time) to bright yellow/green (later in time). The range between the earliest and latest time-delay values are significantly shorter for lag projections on global-signal regressed data.
# %% [markdown]
# # Supplementary Figure G - Consistency in Zero-lag FC Topographies at Finer-Grained Solutions
# %%
base_dir = 'results/cross_val'
fps = [
'caps_smg.dtseries.nii',
'caps_prec.dtseries.nii',
'caps_sm.dtseries.nii',
'eigenmap.dtseries.nii',
'hmm_mean_map.dtseries.nii',
'pca_varimax.dtseries.nii',
'pca.dtseries.nii',
's_ica.dtseries.nii',
't_ica.dtseries.nii'
]
comps_all = []
for val in range(12):
comp_dir=f'comp{val+1}'
comps_val = []
for fp in fps:
_, cifti_maps, _ = pull_cifti_data(load_cifti(f'{base_dir}/{comp_dir}/{fp}'))
comps_val.append(cifti_maps)
comps_all.append(np.vstack(comps_val))
mean_abs_corr = []
for comp_val in comps_all:
corr_mat = np.corrcoef(comp_val)
corr_ltr = corr_mat[np.tril_indices(corr_mat.shape[0],k=1)]
mean_abs_corr.append(np.mean(np.abs(corr_ltr)))
fig, ax = plt.subplots(figsize=(7,7))
ax.set_title('Mean Abs. Correlation by # of Dimensions', fontweight='bold', fontsize=17)
ax.set_xlabel('# of Dimensions', fontsize=13)
ax.set_ylabel('Mean Abs. Correlation', fontsize=13)
ax.plot(range(1,13), mean_abs_corr)
plt.savefig('results/figures/corr_by_number.eps', bbox_inches='tight')
# %% [markdown]
# ## Appendix I - Code to Calculate Average Duration of Complex Principal Components
# %% [markdown]
# The following code was used to calculate the duration of the first three complex principal components. The procedure was as follows: 1) the temporal phase was derived from the complex principal component time series, 2) the phase was unwrapped, and 3) the average duration was calculated as the average time from start (0) to end (2pi) for all cycles of the component time series.
# %%
# pca_complex_res = pickle.load(open('results/pca_rest_complex_results.pkl', 'rb'))
# comp0_phase = np.unwrap(np.angle(pca_complex_res['pc_scores'][:,0]))
# comp1_phase = np.unwrap(np.angle(pca_complex_res['pc_scores'][:,1]))
# comp2_phase = np.unwrap(np.angle(pca_complex_res['pc_scores'][:,2]))
# avg_cycle_comp0 = (comp0_phase[-1]-comp0_phase[0])/60000
# avg_cycle_comp0 = (2*np.pi)/avg_cycle_comp0
# avg_cycle_comp1 = (comp1_phase[-1]-comp1_phase[0])/60000
# avg_cycle_comp1 = (2*np.pi)/avg_cycle_comp1
# avg_cycle_comp2 = (comp2_phase[-1]-comp2_phase[0])/60000
# avg_cycle_comp2 = (2*np.pi)/avg_cycle_comp2
# %% [markdown]
# ## Scratch Code
# %%
pca_complex_res = pickle.load(open('results/pca_rest_complex_results.pkl', 'rb'))
pca_ts = pickle.load(open('demo_files/pca_ts.pkl', 'rb'))
sica_ts = pickle.load(open('demo_files/sica_ts.pkl', 'rb'))
pca_ts_c = pickle.load(open('demo_files/pca_gs_ts.pkl', 'rb'))
varimax_ts = pickle.load(open('demo_files/varimax_ts.pkl', 'rb'))
gs_signal = pickle.load(open('demo_files/gs_results.pkl', 'rb'))
qpp_ts = pickle.load(open('results/qpp_results.pkl', 'rb'))[3]
pca_complex_ts = pca_complex_res['pca']['pc_scores']
xcorr(zscore(np.imag(pca_complex_ts[:,0])), zscore(np.real(pca_complex_ts[:,2]).T), maxlags=30) #
# %%
_, cifti_maps_real, n_time = pull_cifti_data(load_cifti('results/pca_rest_complex_real.dtseries.nii'))
_, cifti_maps_imag, n_time = pull_cifti_data(load_cifti('results/pca_rest_complex_imag.dtseries.nii'))
_, cifti_eigenmap_p90, n_time = pull_cifti_data(load_cifti('results/eigenmap_thres/eigenmap_p90.dtseries.nii'))
_, cifti_maps_pseed, n_time = pull_cifti_data(load_cifti('results/fc_map_precuneus.dtseries.nii'))
_, cifti_maps_pseed_gs, n_time = pull_cifti_data(load_cifti('results/fc_map_gs_precuneus.dtseries.nii'))
_, cifti_maps_smseed, n_time = pull_cifti_data(load_cifti('results/fc_map_sm.dtseries.nii'))
_, cifti_maps_smseed_gs, n_time = pull_cifti_data(load_cifti('results/fc_map_gs_sm.dtseries.nii'))
_, cifti_maps_spseed, n_time = pull_cifti_data(load_cifti('results/fc_map_supramarginal.dtseries.nii'))
_, cifti_maps_spseed_gs, n_time = pull_cifti_data(load_cifti('results/fc_map_gs_supramarginal.dtseries.nii'))
zero_mask = np.std(cifti_maps_real, axis=0) > 0
cifti_maps_real = cifti_maps_real[:, zero_mask].copy()
cifti_maps_imag = cifti_maps_imag[:, zero_mask].copy()
cifti_eigenmap_p90 = cifti_eigenmap_p90[0, zero_mask].copy()
cifti_maps_pseed = cifti_maps_pseed[:, zero_mask].copy()
cifti_maps_pseed_gs = cifti_maps_pseed_gs[:, zero_mask].copy()
cifti_maps_smseed = cifti_maps_smseed[:, zero_mask].copy()
cifti_maps_smseed_gs = cifti_maps_smseed_gs[:, zero_mask].copy()
cifti_maps_spseed = cifti_maps_spseed[:, zero_mask].copy()
cifti_maps_spseed_gs = cifti_maps_spseed_gs[:, zero_mask].copy()

figures.ipynb at commit 96e91dd, no license · at the source

Overview

Authors: Shuo Lv¹, Jinlong Li¹, Ruoqi Yang¹, Xinyu Wu¹, Zhiming Wang¹, Wenjing Zhu¹, Tan Gao¹, Jia-Hong Gao², Guoyuan Yang¹

ORCID iDs: Shuo Lv, Xinyu Wu, Jia-Hong Gao, Guoyuan Yang

Institutions: Beijing Institute of Technology (China); Peking University (China)

Journal: Nature communications, volume 17, issue 1, article 6232

Dates: received 7 September 2025; accepted 24 April 2026; published online 8 May 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41467-026-72931-6 · PMID 42103777 · PMCID PMC13369909 · OpenAlex W7160623821

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: fMRI (modality), human (organism), cognitive (subfield)

Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Graphs, fMRI & imaging

Keywords: Dynamical systems, Cognitive neuroscience

MeSH: Cerebellum*, Adult, Connectome, Female, Humans, Magnetic Resonance Imaging, Male, Nerve Net, Neural Pathways, Principal Component Analysis (* major topic)

Topic: Vestibular and auditory disorders (Neurology, Neuroscience), according to OpenAlex

Funding: National Natural Science Foundation of China (National Science Foundation of China) (82302175, 62336002)

Citations: not cited yet (Europe PMC); 120 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 18 matches between paragraphs and lines of code.

RaichleLab/lag-code

License: other

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: e8b6e17d42d1808476bdcdc81230291b698a4f6c, 6 November 2024

Languages: MATLAB (8)

Size: 10 files, 8 scripts

Software Heritage: not archived

Found in: the text, “Lag projection”

Holds: README, license file

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Statistics and Machine Learning Toolbox (3 files), Signal Processing Toolbox (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

10 files

acf_hwhm.m, MATLAB, 35 lines
create_blocks.m, MATLAB, 20 lines
f_alpha_gaussian.m, MATLAB, 107 lines
lagged_cov.m, MATLAB, 30 lines
parabolic_interp.m, MATLAB, 49 lines
spectral_template.m, MATLAB, 88 lines
surrogate_TDE.m, MATLAB, 167 lines
tdmx_template.m, MATLAB, 153 lines
LICENSE.txt, License, 21 lines
README.md, Text, 26 lines

At the source: github.com/RaichleLab/lag-code

BIT-YangLab/CPCA_Cerebellum

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: db528107547c07f91093022b5b46532c4a7b2d09, 2 April 2026

Languages: Python (35), Shell (1)

Size: 446 files, 36 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (34 files), SciPy (31 files), NiBabel (25 files), Matplotlib (20 files), pandas (7 files), scikit-learn (7 files), BrainSpace (2 files), seaborn (2 files), statsmodels (2 files), PyTorch (1 file), Connectome Workbench (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

37 files

BOLD_simulation.py, Python, 412 lines, 1 match
CPCA_analysis.py, Python, 508 lines, 2 matches
ROC_CM_null.py, Python, 158 lines, 2 matches
cpca_cere.py, Python, 153 lines
cpca_cere_subjs.py, Python, 211 lines
cpca_cortex.py, Python, 193 lines
cpca_cortex_subjs.py, Python, 211 lines
cpca_full.py, Python, 197 lines
cpca_full_subjs.py, Python, 201 lines
exp_var.py, Python, 46 lines
hmm_cere.py, Python, 97 lines
moran.py, Python, 169 lines
phase_crp.py, Python, 69 lines
phase_crp_cortex.py, Python, 69 lines
phase_crp_full.py, Python, 65 lines
run_fc_matrix.py, Python, 94 lines
sCCA_analysis_wholebrain.py, Python, 716 lines, 4 matches
sex_LSVM.py, Python, 463 lines, 2 matches
sex_volume.py, Python, 105 lines
sex_volume_hist.py, Python, 227 lines
sex_wholebrain_hist.py, Python, 381 lines
utils/__init__.py, Python, 1 line
utils/complex_fastica.py, Python, 181 lines
utils/giftis_to_cifti.sh, Shell, 11 lines
utils/rotation.py, Python, 135 lines
utils/utils.py, Python, 144 lines
utils/utils_cere.py, Python, 103 lines
utils/utils_cere_subjs.py, Python, 135 lines
utils/utils_cortex.py, Python, 102 lines
utils/utils_cortex_subjs.py, Python, 135 lines
utils/utils_cortex_subjs_srp.py, Python, 135 lines
utils/utils_full.py, Python, 104 lines
utils/utils_full_down_add.py, Python, 168 lines
utils/utils_full_down_subjs.py, Python, 168 lines
utils/utils_full_subjs.py, Python, 135 lines
utils/utils_full_subjs_srp.py, Python, 135 lines
README.md, Text, 1 line

At the source: github.com/BIT-YangLab/CPCA_Cerebellum

tsb46/BOLD_WAVES

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 96e91ddccdacee4f6a0faa623d28528937fc9313, 5 September 2023

Languages: Python (17), Shell (9), Jupyter (2), R (1)

Size: 202 files, 29 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, environment (requirements.txt, requirements_notebook.txt), 2 notebooks

Not found: license file, CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (19 files), SciPy (16 files), scikit-learn (9 files), NiBabel (5 files), Connectome Workbench (5 files), Matplotlib (3 files), pandas (2 files), Brain Connectivity Toolbox (1 file), BrainSpace (1 file), Nilearn (1 file), statsmodels (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

30 files

data/get_rest_cifti.R, R, 23 lines
figures.ipynb, Jupyter, 2,209 lines, 5 matches
preprocess/fsaverage_resample.sh, Shell, 37 lines
preprocess/global_signal_regress.py, Python, 82 lines
preprocess/global_signal_regress.sh, Shell, 14 lines
preprocess/merge_metric_resample.sh, Shell, 35 lines
preprocess/norm_filter.py, Python, 74 lines
preprocess/norm_filter.sh, Shell, 8 lines
preprocess/parcel_ts.sh, Shell, 12 lines
preprocess/surface_smooth.sh, Shell, 11 lines
preprocess_rest.sh, Shell, 38 lines
run_analysis.sh, Shell, 95 lines, 1 match
run_cap_analysis.py, Python, 187 lines
run_cpca_reconstruction.py, Python, 70 lines
run_eigenmap.py, Python, 111 lines
run_fc_matrix.py, Python, 83 lines
run_global_signal_analysis.py, Python, 74 lines
run_hmm.py, Python, 87 lines
run_ica.py, Python, 119 lines
run_lag_projection.py, Python, 169 lines
run_main_pca.py, Python, 194 lines
run_peak_average.py, Python, 133 lines
run_qpp.py, Python, 262 lines, 1 match
run_seed_fc.py, Python, 110 lines
simulation.ipynb, Jupyter, 898 lines
utils/complex_fastica.py, Python, 181 lines
utils/giftis_to_cifti.sh, Shell, 11 lines
utils/rotation.py, Python, 135 lines
utils/utils.py, Python, 153 lines
README.md, Text, 126 lines

At the source: github.com/tsb46/BOLD_WAVES

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: BIT-YangLab/CPCA_Cerebellum, tsb46/BOLD_WAVES

Read it in the paper: doi.org/10.1038/s41467-026-72931-6.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
73 scripts, each with its path and the digest of its content;
18 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

github.com/esfinn/cpm_tutorial, at github.com; found in the text, “Behavioral phenotypes are associated with…”
humanconnectome.org/study/hcp-young-adult, at Human Connectome Project; found in “Data availability”

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: humanconnectome.org/study/hcp-young-adult
it says that the data are available on request

Read it in the paper: doi.org/10.1038/s41467-026-72931-6.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 2 keywords, 10 MeSH terms, 1 funder, 118 references.

Cite

This paper

Lv, S., Li, J., Yang, R., Wu, X., Wang, Z., Zhu, W., Gao, T., Gao, J.-H., & Yang, G. (2026). Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior. Nature communications, 17(1), 6232. https://doi.org/10.1038/s41467-026-72931-6

BibTeX

@article{lv2026three, author = {Lv, Shuo and Li, Jinlong and Yang, Ruoqi and Wu, Xinyu and Wang, Zhiming and Zhu, Wenjing and Gao, Tan and Gao, Jia-Hong and Yang, Guoyuan}, title = {{Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior}}, journal = {Nature communications}, year = {2026}, month = may, volume = {17}, number = {1}, pages = {6232}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-72931-6}, url = {https://doi.org/10.1038/s41467-026-72931-6}, pmid = {42103777}, pmcid = {PMC13369909} }

RIS

TY - JOUR AU - Lv, Shuo AU - Li, Jinlong AU - Yang, Ruoqi AU - Wu, Xinyu AU - Wang, Zhiming AU - Zhu, Wenjing AU - Gao, Tan AU - Gao, Jia-Hong AU - Yang, Guoyuan TI - Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/05/08 VL - 17 IS - 1 SP - 6232 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-72931-6 UR - https://doi.org/10.1038/s41467-026-72931-6 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-72931-6", "type": "article-journal", "title": "Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior", "container-title": "Nature communications", "author": [ { "family": "Lv", "given": "Shuo" }, { "family": "Li", "given": "Jinlong" }, { "family": "Yang", "given": "Ruoqi" }, { "family": "Wu", "given": "Xinyu" }, { "family": "Wang", "given": "Zhiming" }, { "family": "Zhu", "given": "Wenjing" }, { "family": "Gao", "given": "Tan" }, { "family": "Gao", "given": "Jia-Hong" }, { "family": "Yang", "given": "Guoyuan" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "6232", "DOI": "10.1038/s41467-026-72931-6", "PMID": "42103777", "PMCID": "PMC13369909", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-72931-6", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 8 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76011-7 [code]: Human cortex organizes dynamic co-fluctuations along the sensorimotor-association axis.
Journal: Nature communications
In common: BrainSpace, Brain Connectivity Toolbox, Connectome Workbench, 7 other tools, 13 references
[2] doi:10.21203/rs.3.rs-9326213/v1 [code]: Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: Connectome Workbench, Nilearn, NiBabel, 7 other tools, humanconnectome.org/study/hcp-young-adult, fMRI, cognitive, 11 references
[3] doi:10.64898/2026.03.09.710558 [code]: Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: bioRxiv (preprint)
In common: Connectome Workbench, Nilearn, NiBabel, 7 other tools, humanconnectome.org/study/hcp-young-adult, fMRI, cognitive, 11 references
[4] doi:10.1038/s41467-026-72940-5 [code]: Cerebellar growth is associated with domain-specific cerebral maturation and socio-linguistic behavior.
Journal: Nature communications
In common: Connectome Workbench, NiBabel, statsmodels, 6 other tools, 13 references
[5] doi:10.1002/hbm.70483 [code]: Untamed: Unconstrained Tensor Decomposition and Graph Node Embedding for Cortical Parcellation.
Journal: Human brain mapping
In common: Connectome Workbench, Nilearn, Signal Processing Toolbox, 8 other tools, humanconnectome.org/study/hcp-young-adult, fMRI, 4 references
[6] doi:10.1038/s41467-026-71270-w [code]: Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan.
Journal: Nature communications
In common: BrainSpace, Connectome Workbench, Nilearn, 8 other tools, fMRI, 3 references, author Jia-Hong Gao
[7] doi:10.1016/j.isci.2026.116903 [code]: Neurobiological and behavioral relevance of intrinsic functional connectome constraints on task-evoked neural activation.
Journal: iScience
In common: Connectome Workbench, Nilearn, NiBabel, 6 other tools, humanconnectome.org/study/hcp-young-adult, fMRI, cognitive, 4 references
[8] doi:10.1038/s41467-026-71151-2 [code]: Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: Brain Connectivity Toolbox, Connectome Workbench, Nilearn, 10 other tools, cognitive, 3 references
[9] doi:10.1038/s41398-026-04025-2 [code]: Brain energetic landscapes shape state dysregulation in major depressive disorder: a morphological network controllability perspective.
Journal: Translational psychiatry
In common: BrainSpace, Brain Connectivity Toolbox, Connectome Workbench, 10 other tools, 3 references
[10] doi:10.1038/s41467-026-75959-w [code]: Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: BrainSpace, Nilearn, NiBabel, 8 other tools, 6 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

RaichleLab/lag-code

BIT-YangLab/CPCA_Cerebellum

tsb46/BOLD_WAVES

Code availability statement

Tracing map

Data

Datasets cited

Data availability statement

Versions

Version 1, 28 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity