Higher-order connectomics of human brain function reveals local topological signatures of task decoding, individual identification, and behavior

The study examines whether higher-order interactions (HOIs)—group relationships among three or more brain regions—provide advantages over traditional pairwise functional connectivity (FC) in analyzing human fMRI data. Motivated by the limitations of pairwise network models and increasing evidence for HOIs across micro- and macro-scales, the work aims to test a topological, time-resolved method for reconstructing HOI structures and to compare its performance versus classical node- and edge-centric approaches in three applications: task decoding, functional fingerprinting, and brain–behavior associations. The overarching hypothesis is that HOIs, especially when assessed locally, capture informative and behaviorally relevant aspects of brain dynamics that are missed by pairwise models.

Network models of the connectome have traditionally centered on pairwise relationships (FC) between brain regions. However, theoretical and empirical work indicates that HOIs can qualitatively alter system dynamics and are critical for characterizing brain function. At micro-scales, simultaneous group neuronal firing has been observed in animal studies; at macro-scales, HOIs have been inferred from human neuroimaging using information-theoretic (e.g., multivariate mutual information, integrated information decomposition, partial entropy decomposition) and topological data analysis frameworks. Edge-centric approaches extend beyond node-centric FC and have improved identification of overlapping communities and dynamic connectivity, yet the benefits of HOIs for fMRI remain underexplored. Initial fMRI applications suggest higher-order dependencies may encode biomarkers relevant to consciousness states and aging. Recent work has shown temporal HOI features can improve classification in other domains compared to lower-order, edge-centric measures, motivating a systematic comparison in human fMRI.

Dataset: Human Connectome Project (HCP) Release Q3; 100 unrelated adult subjects (54 females, 46 males; mean age 29.1 ± 3.7). Ethics and informed consent in accordance with HCP protocols. Preprocessing: Minimally preprocessed HCP data; GLM regression including detrending (quadratic trends), motion regressors and derivatives, white matter and CSF signals and derivatives, global signal regression and derivative; bandpass filter 0.01–0.15 Hz; voxel-wise fMRI signals averaged to atlas parcels and z-scored. Parcellation: Schaefer cortical atlas (100 regions) plus 16 subcortical and 3 cerebellar regions from HCP (Atlas_ROI2.nii.gz), total N = 119 nodes. Classical measures: BOLD time series (node signals); edge time series computed as element-wise products of pairs of z-scored node signals; FC (Pearson correlations between node pairs); eFC (correlations between edge time series). Higher-order approach: Generalization to k-order co-fluctuations (up to k = 2; triangles). For each set of k+1 z-scored node signals, compute standardized element-wise product over time; apply sign remapping so perfectly concordant group interactions (all positive or all negative) are positive, discordant are negative; define weights w for edges and triangles at each timepoint. Build an instantaneous weighted simplicial complex K^t with edges and triangles; analyze via computational topology. Topological indicators: Global—(i) Hyper-coherence: fraction of violating coherent triangles whose standardized triangle weight exceeds those of constituent edges (identified during filtration); (ii) Hyper-complexity: sliced Wasserstein distance between the persistence diagram of H1 (1D cycles, birth–death pairs along filtration) and the empty diagram, quantifying topological complexity of coherent/decoherent landscapes. Local—(i) Violating triangles list Δ and weights; (ii) Homological scaffold: weighted graph where edge weights equal the frequency of participation in 1D cycles, highlighting links bridging coherent/decoherent structures. Task decoding pipeline: Concatenate first 300 volumes of resting-state fMRI with seven HCP task runs (excluding rest blocks) to form a unified time series. For each local signal type (BOLD, edges, triangles, scaffold), compute time–time recurrence plots (Pearson correlation between timepoints), binarize at the 95th percentile threshold, and apply Louvain community detection. Evaluate task/rest block decoding by element-centric similarity (ECS) between detected communities and ground-truth blocks. Results averaged over 100 subjects (mean across LR and RL encodings); additional thresholds reported in SI. Fingerprinting: Use two resting-state sessions (REST_1_LR and REST_2_LR for test/retest). From temporal signals, derive static connectivity: node FC, eFC, average violating triangle weights, average scaffold weights. Compute an identifiability matrix across subjects (S=100) using Pearson similarity of test/retest connectivity; differential identifiability (Idiff) is the difference between within-subject (diagonal) and between-subject (off-diagonal) similarities. Compare whole-brain vs local (connections involving at least one node in each of seven functional networks: VIS, SM, DA, VA, L, FP, DMN; plus SC). Brain–behavior association: Partial Least Square Correlation (PLSC) between whole-brain or network-specific connectivity features (FC, eFC, triangles, scaffold) and 10 HCP cognitive domains (episodic memory, executive functions, fluid intelligence, language, processing speed, self-regulation/impulsivity, spatial orientation, sustained attention, verbal episodic memory, working memory). For subdomains with multiple raw scores, reduce via PCA to first principal component. Assess significance via permutation testing (1000 permutations, p < 0.05); reliability via bootstrap (1000 resamples, absolute standard score > 2). Quantify covariance explained by summing squared singular values of significant components normalized by total. Robustness: bootstrap 80/100 subjects repeated 100 times; also compare with Canonical Correlation Analysis (reported in SI).

• Global indicators: Hyper-coherence and hyper-complexity showed remarkably similar values across tasks and rest; no significant differences (all pairwise t-tests p > 0.1), indicating limited utility for task decoding at the whole-brain level. • Task decoding (local): Recurrence plots and community detection demonstrated improved temporal discrimination of task and rest blocks with higher-order local indicators versus lower-order ones. Reported ECS values averaged across 100 subjects: BOLD ECS = 0.46; edge ECS = 0.73; triangle ECS = 0.76; scaffold ECS = 0.88, showing progressive improvement with higher-order methods. • Fingerprinting: Whole-brain identifiability scores were similar across methods. Restricting to local connections within functional networks, the triangle method consistently outperformed FC and eFC across nearly all seven networks; scaffold identifiability remained below 9% for all networks (SI Fig. S3). Triangle-based nodal strength showed high inter-subject variability, particularly in interactions between unimodal (visual, somatosensory) and transmodal (DMN, frontoparietal) systems. • Brain–behavior association: At the whole-brain level, all methods explained similar amounts of covariance (~10–20%). At the local network level, higher-order methods (triangles, scaffold) substantially increased explained covariance compared to FC/eFC, with triangles reaching ~80% covariance in somatosensory areas. Cognitive saliences varied across methods and networks, indicating distinct cognitive dimensions captured by higher-order versus pairwise measures. • Robustness: Bootstrap subsampling and alternative multivariate approaches (CCA) corroborated that local higher-order methods outperform lower-order approaches across most functional networks. • Overall: Local higher-order features capture nuanced, spatially specific brain dynamics that enhance task decoding, individual identification, and brain–behavior associations beyond traditional pairwise models.

The findings address the central question of whether inferred higher-order interactions improve fMRI analyses over traditional pairwise connectivity. Globally, higher-order metrics did not confer advantages, suggesting that whole-brain summaries may mask informative localized structure. Locally, however, higher-order indicators—violating triangles and homological scaffolds—provided markedly better temporal decoding of tasks, stronger subject identifiability within functional subsystems, and more robust multivariate associations with behavior. These results imply that the brain’s functional coordination involves spatially specific, synergistic group interactions not captured by pairwise measures, aligning with theories distinguishing redundancy-dominated pairwise patterns from synergistic higher-order contributions. The approach offers a time-resolved window into recruitment of functional groups and dynamic reconfiguration of subsystems across brain states, enhancing understanding of how unimodal sensory and transmodal association regions interact to shape individual fingerprints and behavioral variability.

A topological, time-resolved framework for reconstructing higher-order interactions in fMRI reveals that local group-level features outperform traditional pairwise measures in three key applications: task decoding, functional brain fingerprinting, and brain–behavior association. While global indicators offer limited discrimination, local higher-order structures capture subtle, spatially specific dynamics central to cognition and behavior. The work demonstrates practical gains of higher-order connectomics for human fMRI and motivates methodological advances that embrace the computational challenges of higher-order modeling. Future research directions include reducing computational burdens by reconstructing targeted subsets of k-order interactions, and integrating causal and geometric tools (e.g., Takens’ embeddings, spectral decomposition) to further elucidate temporal dependencies and mechanisms in brain network dynamics.

The approach is computationally intensive: reconstructing co-fluctuation patterns up to order k scales as O(N^k), requiring approximately 5 minutes per HCP subject on a high-end workstation. Analyses were limited to k = 2 (triangles), and direct measurements of true group interactions are not available in human fMRI, necessitating inference from noisy signals. Global higher-order indicators did not outperform pairwise methods, indicating that utility is primarily local. Future work should focus on selecting informative subsets of higher-order interactions to reduce complexity and exploring complementary methods to probe causality and temporal dependencies.