combinatorial-complex-brain-fmri

name	combinatorial-complex-brain-fmri
description	The Human Brain as a Combinatorial Complex - framework for constructing combinatorial complexes from fMRI time series that captures both pairwise and higher-order neural interactions through information-theoretic measures. Bridges topological deep learning and network neuroscience.

The Human Brain as a Combinatorial Complex

Framework for constructing combinatorial complexes (CCs) from fMRI time series data that captures both pairwise and higher-order neural interactions, bridging topological deep learning and network neuroscience.

Metadata

• Source: arXiv:2511.20692v2
• Authors: Valentina Sánchez, Çiçek Güven, Koen Haak, Theodore Papamarkou, Gonzalo Nápoles
• Published: 2025-11-22
• Category: q-bio.NC, cs.LG, math.AT

Core Methodology

Key Innovation

Traditional graph-based representations of brain networks systematically miss higher-order dependencies that characterize neural complexity. This methodology introduces combinatorial complexes (CCs) as a unified framework that captures:

• Pairwise interactions (edges)
• Higher-order interactions (simplices, hyperedges)
• Information-theoretic relationships between brain regions

Theoretical Foundation

From Graphs to Combinatorial Complexes

Graph Limitations:

• Edges only capture pairwise relationships: G = (V, E)
• Cannot represent simultaneous activity of multiple brain regions
• Misses synergistic interactions (e.g., three regions co-activating)

Combinatorial Complex Solution:

• 0-cells: Vertices (individual brain regions)
• 1-cells: Edges (pairwise interactions)
• k-cells: k-simplices (k+1-way interactions)
• Incidence relations: How cells connect across dimensions

Information-Theoretic Construction

Step 1: Time Series to Information Measures

For each subset S ⊆ V of brain regions:
    Compute mutual information I(S) = H(S) - Σ H(v) for v ∈ S
    or multi-information (total correlation)

Step 2: Threshold-Based Cell Construction

For k = 1 to max_order:
    For each (k+1)-subset S of regions:
        If I(S) > θ_k:
            Add k-cell with vertices S to complex

Step 3: Incidence Structure

• Define boundary operator ∂k: C_k → C{k-1}
• ∂_k maps k-cells to their (k-1)-dimensional faces
• Forms chain complex: ... → C_2 → C_1 → C_0 → 0

Combinatorial Complex Types

1. Simplicial Complex (SC)

• Property: Closed under taking subsets
• Interpretation: If regions {A,B,C} co-activate, all pairs {A,B}, {B,C}, {A,C} also interact
• Advantage: Mathematically well-understood, persistent homology available
• Constraint: May overcount (strong closure property)

2. Cell Complex (CC)

• Property: More flexible incidence relations
• Interpretation: Allows non-simplicial shapes (e.g., cycles without filling)
• Advantage: Better fits neural topology
• Challenge: More complex computational structure

3. Hypergraph

• Property: No closure requirement
• Interpretation: Only specific higher-order interactions captured
• Advantage: Most flexible
• Trade-off: Less algebraic structure for analysis

Deep Learning Integration

Combinatorial Complex Neural Networks (CCNN)

Message Passing on CCs:

For each cell c in dimension k:
    Aggregate messages from boundary ∂(c) and coboundary δ(c)
    Update representation: h_c^{(t+1)} = UPDATE(h_c^{(t)}, AGG({h_b for b ∈ N(c)}))

Higher-Order Convolutions:

• 0-order: Standard node-level features
• 1-order: Edge-level interactions
• k-order: k-way synergistic patterns

Brain Network Applications

fMRI Pipeline:

Preprocessed fMRI → ROI Time Series 
    → Information Estimation 
    → CC Construction 
    → CCNN Learning 
    → Clinical Prediction

Information Estimators:

• Gaussian: For linear correlations
• K-nearest neighbors: Non-parametric mutual information
• Kernel density: Smooth probability estimates

Implementation Guide

Prerequisites

# Required libraries
pip install gudhi            # Topological data analysis
pip install torch-geometric   # Graph neural networks
pip install dit              # Information theory (discrete)
pip install nilearn          # fMRI processing
pip install scipy

Step-by-Step Implementation

Step 1: fMRI Data Preparation

import numpy as np
from nilearn import datasets
from nilearn.connectome import ConnectivityMeasure

def load_fmri_data(n_subjects=100):
    """Load and preprocess resting-state fMRI data."""
    # Use publicly available dataset
    dataset = datasets.fetch_abide_pcp(
        data_dir='./data',
        n_subjects=n_subjects,
        pipeline='cpac',
        band_pass_filter=True
    )
    
    # Extract time series using Schaefer atlas
    atlas = datasets.fetch_atlas_schaefer_2018(n_rois=200)
    
    time_series = []
    for func_file in dataset.func_preproc:
        from nilearn.maskers import NiftiLabelsMasker
        masker = NiftiLabelsMasker(
            labels_img=atlas.maps,
            standardize=True,
            detrend=True
        )
        ts = masker.fit_transform(func_file)
        time_series.append(ts)
    
    return time_series, atlas

Step 2: Information-Theoretic Cell Construction

from sklearn.neighbors import NearestNeighbors
from scipy.stats import entropy

def estimate_mutual_information_knn(X, Y, k=5):
    """K-NN based mutual information estimation."""
    # Kraskov-Stögbauer-Grassberger estimator
    n = len(X)
    
    # Joint space
    XY = np.column_stack([X, Y])
    
    # Find k-nearest neighbors
    nbrs = NearestNeighbors(n_neighbors=k+1).fit(XY)
    distances, _ = nbrs.kneighbors(XY)
    epsilon = distances[:, k]  # Distance to k-th neighbor
    
    # Count neighbors in marginal spaces
    nx = np.array([np.sum(np.abs(X - X[i]) < epsilon[i]) - 1 for i in range(n)])
    ny = np.array([np.sum(np.abs(Y - Y[i]) < epsilon[i]) - 1 for i in range(n)])
    
    # MI estimate
    mi = np.mean(np.log(n / (nx * ny))) + np.log(k) + np.euler_gamma
    
    return max(0, mi)

def construct_simplicial_complex(time_series, threshold_percentile=90, max_dim=3):
    """
    Build simplicial complex from fMRI time series.
    
    Args:
        time_series: [n_regions, n_timepoints] array
        threshold_percentile: Percentile for edge threshold
        max_dim: Maximum simplex dimension
    
    Returns:
        simplices: List of (dimension, vertices) tuples
    """
    n_regions = time_series.shape[0]
    
    # Compute pairwise mutual information
    mi_matrix = np.zeros((n_regions, n_regions))
     i  (n_regions):
         j  (i+, n_regions):
            mi = estimate_mutual_information_knn(
                time_series[i], time_series[j]
            )
            mi_matrix[i, j] = mi_matrix[j, i] = mi
    
    
    threshold = np.percentile(mi_matrix[mi_matrix > ], threshold_percentile)
    
    
    edges = []
     i  (n_regions):
         j  (i+, n_regions):
             mi_matrix[i, j] > threshold:
                edges.append((i, j))
    
    
    simplices = [(, (i,))  i  (n_regions)]  
    simplices.extend([(, edge)  edge  edges])   
    
    
     itertools  combinations
    
     k  (, max_dim + ):
        
         vertices  combinations((n_regions), k+):
            
            is_clique = (
                ((v1, v2), (v1, v2))  edges
                 v1, v2  combinations(vertices, )
            )
             is_clique:
                
                sub_matrix = mi_matrix[np.ix_(vertices, vertices)]
                avg_mi = np.mean(sub_matrix[np.triu_indices_from(sub_matrix, k=)])
                 avg_mi > threshold * :  
                    simplices.append((k, vertices))
    
     simplices, mi_matrix

Step 3: Combinatorial Complex Neural Network

import torch
import torch.nn as nn
import torch.nn.functional as F

class CCNNLayer(nn.Module):
    """Combinatorial Complex Neural Network Layer."""
    
    def __init__(self, in_channels, out_channels, max_dim=3):
        super().__init__()
        self.max_dim = max_dim
        
        # Separate convolutions for each dimension
        self.convs = nn.ModuleList([
            nn.Linear(in_channels, out_channels)
            for _ in range(max_dim + 1)
        ])
        
        # Inter-dimensional message passing
        self.boundary_convs = nn.ModuleDict()
        for k in range(1, max_dim + 1):
            self.boundary_convs[str(k)] = nn.Linear(in_channels, out_channels)
    
    def forward(self, x_dict, incidence_dict):
        """
        Args:
            x_dict: Dictionary {dim: features [n_cells, in_channels]}
            incidence_dict: Dictionary of boundary/coboundary operators
        Returns:
            out_dict: Updated features
        """
        out_dict = {}
        
        for k in range(self.max_dim + 1):
            h = self.convs[k](x_dict[k])
            
            # Aggregate from boundary (lower dimension)
             k >   (k)  incidence_dict:
                boundary_msg = .aggregate_from_boundary(
                    x_dict[k-], incidence_dict[(k)]
                )
                h = h + .boundary_convs[(k)](boundary_msg)
            
            
             k < .max_dim  (k+)  incidence_dict:
                coboundary_msg = .aggregate_from_coboundary(
                    x_dict[k+], incidence_dict[(k+)]
                )
                h = h + coboundary_msg
            
            out_dict[k] = F.relu(h)
        
         out_dict
    
     ():
        
        
         torch.matmul(incidence_matrix, x_lower)
    
     ():
        
        
         torch.matmul(incidence_matrix.T, x_higher)


 (nn.Module):
    
    
     ():
        ().__init__()
        .max_dim = max_dim
        
        
        .encoders = nn.ModuleList([
            nn.Sequential(
                CCNNLayer(feature_dims[i], feature_dims[i+], max_dim),
                CCNNLayer(feature_dims[i+], feature_dims[i+], max_dim)
            )
             i  ((feature_dims) - )
        ])
        
        
        .classifier = nn.Sequential(
            nn.Linear(feature_dims[-] * (max_dim + ), ),
            nn.ReLU(),
            nn.Dropout(),
            nn.Linear(, num_classes)
        )
    
     ():
        
         encoder  .encoders:
            x_dict = encoder(x_dict, incidence_dict)
        
        
        pooled = []
         k  (.max_dim + ):
             k  x_dict:
                pooled.append(x_dict[k].mean(dim=))
        
        
        x = torch.cat(pooled, dim=-)
        
         .classifier(x)

Step 4: Persistent Homology (Optional Enhancement)

import gudhi

def compute_persistent_features(simplices, max_dim=3):
    """Compute persistent homology features from simplicial complex."""
    # Build simplex tree
    st = gudhi.SimplexTree()
    
    for dim, vertices in simplices:
        st.insert(list(vertices), filtration=0.0)
    
    # Compute persistence
    st.compute_persistence()
    
    # Extract features
    features = {}
    for dim in range(max_dim + 1):
        persistence = st.persistence_intervals_in_dimension(dim)
        if len(persistence) > 0:
            features[f'pers_entropy_dim{dim}'] = -np.sum(
                (persistence[:, 1] - persistence[:, 0]) / 
                np.sum(persistence[:, 1] - persistence[:, 0]) * 
                np.log((persistence[:, 1] - persistence[:, 0]) / 
                       np.sum(persistence[:, 1] - persistence[:, 0]))
            )
            features[f'num_features_dim{dim}'] = len(persistence)
    
    return features

Step 5: Training Pipeline

def train_brain_ccnn(model, train_loader, val_loader, epochs=100):
    """Training loop for brain CCNN."""
    optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
    criterion = nn.CrossEntropyLoss()
    
    best_val_acc = 0
    for epoch in range(epochs):
        # Training
        model.train()
        train_loss = 0
        for batch in train_loader:
            x_dict, incidence_dict, labels = batch
            
            optimizer.zero_grad()
            outputs = model(x_dict, incidence_dict)
            loss = criterion(outputs, labels)
            loss.backward()
            optimizer.step()
            
            train_loss += loss.item()
        
        # Validation
        model.eval()
        val_correct = 0
        val_total = 0
        with torch.no_grad():
            for batch in val_loader:
                x_dict, incidence_dict, labels = batch
                outputs = model(x_dict, incidence_dict)
                _, predicted = torch.max(outputs, 1)
                val_total += labels.size(0)
                val_correct += (predicted == labels).sum().item()
        
        val_acc = val_correct / val_total
        if val_acc > best_val_acc:
            best_val_acc = val_acc
            torch.save(model.state_dict(), 'best_brain_ccnn.pt')
        
        if (epoch + 1) % 10 == 0:
            print(f"Epoch {epoch+1}: Loss={train_loss/len(train_loader):f}, Val Acc=")
    
     best_val_acc

Applications

• Brain disorder classification (ADHD, ASD, Alzheimer's)
• Functional connectivity analysis beyond pairwise
• Neural complexity quantification
• Higher-order information flow in brain networks
• Multi-scale brain network dynamics

Advantages Over Graph Methods

Aspect	Graph	Combinatorial Complex
Interactions	Pairwise only	Any order
Structure	Edges	Simplices, cells
Algebra	Adjacency matrix	Boundary operators
Analysis	Spectral	Topological (homology)
Expressiveness	Limited	Rich

Pitfalls

• Computational cost: Higher-order complexes scale combinatorially
• Threshold sensitivity: Cell construction depends on MI threshold
• Interpretation: Higher-order features need neuroscientific validation
• Data requirements: Need sufficient time points for reliable MI estimation
• Dimension choice: Optimal max_dim varies by dataset

Related Skills

• higher-order-brain-networks
• dcho-higher-order-brain-connectivity
• brain-higher-order-structures
• topological-quantum-computing

References

@article{sanchez2025combinatorial,
  title={The Human Brain as a Combinatorial Complex},
  author={Sánchez, Valentina and Güven, Çiçek and Haak, Koen and Papamarkou, Theodore and Nápoles, Gonzalo},
  journal={arXiv preprint arXiv:2511.20692},
  year={2025}
}