By name: geometry-aware-brain-dynamics-mapping description: "Geometry-Aware Framework for noninvasive whole human brain dynamics mapping. Incorporates individual cortical geometry into electrophysiology for accurate spatiotemporal reconstruction of brain activity. Activation: geometry-aware, brain dynamics mapping, cortical geometry, noninvasive electrophysiology."
Geometry-Aware Framework for Noninvasive Whole Human Brain Dynamics Mapping
A geometry-aware framework that incorporates individual cortical geometry into noninvasive electrophysiology to accurately reconstruct whole-brain spatiotemporal dynamics from EEG/MEG.
Metadata
- Source: arXiv:2604.25592v1
- Authors: Youssuf Saleh, Camille Gontier, Jonathan Arreguit, Denis Rivière, Pamela Villagrán, Bertrand Thirion, Alain Destexhe
- Published: 2026-04-28
- Categories: q-bio.NC, cs.LG
Core Methodology
Problem Statement
Non-invasive electrophysiology (EEG/MEG) lacks methods that accurately reconstruct whole-brain spatiotemporal dynamics while incorporating individual cortical geometry. Current approaches suffer from:
- Spatial Smearing: Standard source localization blurs activity across regions
- Geometry Neglect: Flattening or ignoring individual cortical folding patterns
- Temporal Limitations: Poor temporal resolution in fMRI-based methods
- Inverse Problem Ill-Posedness: Non-unique solutions in EEG/MEG source localization
Key Innovation
Geometry-Aware Framework integrates:
- Individual Cortical Geometry: Uses subject-specific cortical surface meshes
- Geometric Deep Learning: GNN/CNN architectures operating on cortical surfaces
- Physics-Informed Regularization: Enforces biophysical constraints from electromagnetic theory
- Spatiotemporal Modeling: Joint spatial and temporal dynamics reconstruction
Technical Framework
1. Cortical Geometry Representation
Cortical Surface Representation:
├── High-resolution triangular mesh (vertex ~1-2mm)
├── Sulcal-gyral pattern encoding
├── Local curvature features
├── Distance along cortical surface (geodesic)
└── Normal vectors for orientation
Individual Geometry Processing:
1. T1-weighted MRI → Freesurfer reconstruction
2. White matter surface extraction
3. Mesh decimation (100k-200k vertices)
4. Geodesic distance matrix computation
5. Curvature and thickness feature extraction
2. Forward Physics Model
The EEG/MEG forward problem maps neural currents to sensor measurements:
Y = L · J + ε
Where:
- Y: Sensor measurements [n_sensors × n_timepoints]
- L: Leadfield matrix [n_sensors × n_sources]
- J: Source currents [n_sources × n_timepoints]
- ε: Noise
Leadfield computation (boundary element method):
L_ij = ∫ G(r_sensor_i, r_source_j) · n_j dA
Where G is the Green's function for the head model
3. Geometry-Aware Inverse Solver
J* = argmin_J ||Y - L·J||² + λ₁·R_geometry(J) + λ₂·R_temporal(J) + λ₃·R_spatial(J)
Geometry Regularization:
R_geometry(J) = Σᵢ Σⱼ wᵢⱼ · ||Jᵢ - Jⱼ||² · d_geodesic(i,j)⁻¹
Where wᵢⱼ encodes cortical connectivity strength
and d_geodesic is geodesic distance along cortex
4. Graph Neural Network Architecture
# Geometry-Aware GNN
Input: Source activity J on cortical mesh vertices
GNN Layers:
1. Spatial Graph Convolution:
J' = σ( D⁻¹ᐟ² A D⁻¹ᐟ² J W )
Where A is adjacency matrix based on geodesic distance
D is degree matrix
W is learnable weights
σ is activation (e.g., ReLU)
2. Temporal Convolution:
J'' = TemporalConv(J', kernel_size=K)
3. Geometry-Aware Pooling:
- Pool across hierarchical cortical parcellations
- Preserve sulcal-gyral boundaries
Output: Reconstructed spatiotemporal dynamics
Implementation Guide
Prerequisites
- Python >= 3.8
- PyTorch >= 1.10
- PyTorch Geometric (for GNNs)
- Freesurfer (for cortical reconstruction)
- MNE-Python (for EEG/MEG processing)
- NumPy, SciPy
Step-by-Step Implementation
1. Cortical Geometry Processing
import nibabel as nib
import numpy as np
from scipy.sparse import csr_matrix
from scipy.spatial.distance import cdist
class CorticalGeometry:
"""
Process individual cortical geometry for brain mapping
"""
def __init__(self, surface_file, curvature_file=None):
"""
Initialize with Freesurfer surface
Args:
surface_file: Path to .pial or .white surface file
curvature_file: Optional curvature map
"""
# Load surface mesh
self.surf = nib.freesurfer.read_geometry(surface_file)
self.vertices = self.surf[0] # [n_vertices, 3]
self.faces = self.surf[1] # [n_faces, 3]
self.n_vertices = len(self.vertices)
# Compute geodesic distances
self.geo_dist = self._compute_geodesic_distance()
# Curvature features
if curvature_file:
self.curvature = nib.freesurfer.read_morph_data(curvature_file)
else:
self.curvature = self._estimate_curvature()
# Surface normals
self.normals = self._compute_normals()
def _compute_geodesic_distance(self, max_dist=50):
"""
Compute approximate geodesic distances on cortical surface
"""
from scipy.sparse.csgraph import dijkstra
# Build adjacency from faces
adj = np.zeros((self.n_vertices, self.n_vertices))
for face in self.faces:
for i in range(3):
for j in range(i+1, 3):
v1, v2 = face[i], face[j]
dist = np.linalg.norm(self.vertices[v1] - self.vertices[v2])
adj[v1, v2] = dist
adj[v2, v1] = dist
# Compute geodesic distances (floyd-warshall or dijkstra)
geo_dist = dijkstra(adj, directed=False, limit=max_dist)
return geo_dist
def _estimate_curvature(self):
"""
Estimate mean curvature at each vertex
"""
curvature = np.zeros(self.n_vertices)
for i in range(self.n_vertices):
# Get neighbors
neighbors = np.unique(self.faces[
np.any(self.faces == i, axis=1)
].flatten())
neighbors = neighbors[neighbors != i]
if len(neighbors) > 0:
# Simple curvature estimate based on local variation
local_pts = self.vertices[neighbors]
centroid = np.mean(local_pts, axis=0)
curvature[i] = np.linalg.norm(self.vertices[i] - centroid)
return curvature
def _compute_normals(self):
"""
Compute vertex normals as average of face normals
"""
normals = np.zeros((self.n_vertices, 3))
for face in self.faces:
v1, v2, v3 = self.vertices[face]
fnormal = np.cross(v2 - v1, v3 - v1)
fnormal = fnormal / (np.linalg.norm(fnormal) + 1e-10)
for v in face:
normals[v] += fnormal
# Normalize
normals = normals / (np.linalg.norm(normals, axis=1, keepdims=True) + 1e-10)
return normals
def get_geometry_features(self):
"""
Return geometry features for each vertex
"""
return {
'vertices': self.vertices,
'curvature': self.curvature,
'normals': self.normals,
'geodesic_distance': self.geo_dist,
'sulcal_depth': self._estimate_sulcal_depth()
}
def _estimate_sulcal_depth(self):
"""
Estimate sulcal depth as proxy for folding
"""
# Simple proxy: distance from convex hull
from scipy.spatial import ConvexHull
hull = ConvexHull(self.vertices)
# For each vertex, distance to nearest hull point
hull_pts = self.vertices[hull.vertices]
depths = np.min(cdist(self.vertices, hull_pts), axis=1)
return depths
2. Leadfield Computation
import mne
import numpy as np
class ForwardModel:
"""
Compute EEG/MEG forward model with geometry
"""
def __init__(self, info, trans, bem, src):
"""
Args:
info: MNE info structure with sensor positions
trans: MRI-to-head transformation
bem: Boundary element model
src: Source space (cortical surface)
"""
self.info = info
self.trans = trans
self.bem = bem
self.src = src
# Compute leadfield
self.leadfield = self._compute_leadfield()
def _compute_leadfield(self):
"""
Compute leadfield matrix using BEM
"""
from mne.forward import make_forward_solution
fwd = make_forward_solution(
self.info, self.trans, self.src, self.bem,
meg=True, eeg=True, mindist=5.0
)
# Extract leadfield as dense matrix
leadfield = fwd['sol']['data'] # [n_sensors, n_sources*3]
# For simplified orientation (normal to cortex)
n_sources = len(self.src[0]['vertno'])
leadfield_normal = np.zeros((leadfield.shape[0], n_sources))
for i in range(n_sources):
# Get orientation (normal to surface)
nn = self.src[0]['nn'][i]
# Extract normal component
idx = slice(i*3, (i+1)*3)
leadfield_normal[:, i] = leadfield[:, idx] @ nn
return leadfield_normal
def project_to_sensors(self, source_activity):
"""
Project source activity to sensor space
Args:
source_activity: [n_sources, n_times]
Returns:
sensor_data: [n_sensors, n_times]
"""
return self.leadfield @ source_activity
def get_snr_weights(self, noise_cov):
"""
Compute SNR-based sensor weighting
"""
noise_var = np.diag(noise_cov)
weights = 1.0 / (noise_var + 1e-10)
return weights
3. Geometry-Aware GNN
import torch
import torch.nn as nn
from torch_geometric.nn import GCNConv, GATConv, global_mean_pool
from torch_geometric.data import Data
class GeometryAwareGNN(nn.Module):
"""
Graph Neural Network operating on cortical surface
"""
def __init__(self, in_channels=1, hidden_channels=64, num_layers=4):
super().__init__()
self.convs = nn.ModuleList()
self.batch_norms = nn.ModuleList()
# First layer
self.convs.append(GATConv(in_channels, hidden_channels, heads=4, concat=False))
self.batch_norms.append(nn.BatchNorm1d(hidden_channels))
# Hidden layers
for _ in range(num_layers - 1):
self.convs.append(GATConv(hidden_channels, hidden_channels, heads=4, concat=False))
self.batch_norms.append(nn.BatchNorm1d(hidden_channels))
self.temporal_conv = nn.Conv1d(hidden_channels, hidden_channels, kernel_size=3, padding=1)
def forward(self, x, edge_index, edge_attr=None, batch=None):
"""
Args:
x: Node features [n_nodes, in_channels, n_time]
edge_index: Graph connectivity [2, n_edges]
edge_attr: Edge weights (geodesic distance) [n_edges]
"""
n_time = x.shape[2]
outputs = []
# Process each timepoint
for t in range(n_time):
x_t = x[:, :, t] # [n_nodes, in_channels]
# Graph convolutions
for conv, bn in zip(self.convs, self.batch_norms):
x_t = conv(x_t, edge_index, edge_attr)
x_t = bn(x_t)
x_t = torch.relu(x_t)
outputs.append(x_t)
# Stack temporal
x = torch.stack(outputs, dim=2) # [n_nodes, hidden_channels, n_time]
# Temporal convolution
x = self.temporal_conv(x)
return x
class BrainDynamicsReconstructor(nn.Module):
"""
End-to-end brain dynamics reconstruction
"""
def __init__(self, n_sensors, n_sources, geometry_info):
super().__init__()
self.n_sensors = n_sensors
self.n_sources = n_sources
# Sensor to source mapping
self.sensor_encoder = nn.Sequential(
nn.Linear(n_sensors, 256),
nn.ReLU(),
nn.Linear(256, 128)
)
# Graph network on cortical surface
self.gnn = GeometryAwareGNN(in_channels=128, hidden_channels=64)
# Source decoder
self.source_decoder = nn.Sequential(
nn.Linear(64, 32),
nn.ReLU(),
nn.Linear(32, 1) # Reconstructed source amplitude
)
# Learned leadfield refinement
self.leadfield_correction = nn.Linear(n_sources, n_sources, bias=False)
# Store geometry
self.edge_index = geometry_info['edge_index']
self.edge_attr = geometry_info.get('edge_attr')
def forward(self, sensor_data):
"""
Reconstruct source activity from sensor data
Args:
sensor_data: [batch, n_sensors, n_time]
Returns:
source_activity: [batch, n_sources, n_time]
"""
batch_size, n_time = sensor_data.shape[0], sensor_data.shape[2]
# Encode sensor data
sensor_encoded = []
for t in range(n_time):
enc = self.sensor_encoder(sensor_data[:, :, t])
sensor_encoded.append(enc)
sensor_encoded = torch.stack(sensor_encoded, dim=2) # [batch, 128, n_time]
# Prepare batch for GNN
all_sources = []
for b in range(batch_size):
# Repeat encoded features for each source
x = sensor_encoded[b].unsqueeze(0).repeat(self.n_sources, 1, 1) # [n_sources, 128, n_time]
# GNN forward
gnn_out = self.gnn(x, self.edge_index, self.edge_attr) # [n_sources, 64, n_time]
# Decode to source activity
source_t = []
for t in range(n_time):
dec = self.source_decoder(gnn_out[:, :, t])
source_t.append(dec)
source_activity = torch.stack(source_t, dim=1) # [n_sources, n_time]
all_sources.append(source_activity)
return torch.stack(all_sources, dim=0) # [batch, n_sources, n_time]
4. Training with Physics Constraints
class GeometryAwareLoss(nn.Module):
"""
Multi-component loss with geometry constraints
"""
def __init__(self, leadfield, geo_dist, lambda_data=1.0,
lambda_smooth=0.1, lambda_geodesic=0.05):
super().__init__()
self.leadfield = leadfield
self.geo_dist = geo_dist
self.lambda_data = lambda_data
self.lambda_smooth = lambda_smooth
self.lambda_geodesic = lambda_geodesic
def forward(self, pred_sources, true_sources, sensor_data):
"""
Compute geometry-aware loss
Args:
pred_sources: [batch, n_sources, n_time]
true_sources: [batch, n_sources, n_time]
sensor_data: [batch, n_sensors, n_time]
"""
# Data fidelity: reconstructed sensors should match
pred_sensors = torch.matmul(self.leadfield, pred_sources)
loss_data = torch.mean((pred_sensors - sensor_data)**2)
# Source prediction loss
loss_source = torch.mean((pred_sources - true_sources)**2)
# Spatial smoothness (geometry-aware)
# Nearby sources on cortex should have similar activity
loss_smooth = self._geometry_smoothness(pred_sources)
# Geodesic distance regularization
loss_geodesic = self._geodesic_regularization(pred_sources)
total_loss = (self.lambda_data * loss_data +
loss_source +
self.lambda_smooth * loss_smooth +
self.lambda_geodesic * loss_geodesic)
return {
'total': total_loss,
'data': loss_data,
'source': loss_source,
'smooth': loss_smooth,
'geodesic': loss_geodesic
}
def _geometry_smoothness(self, sources):
"""
Enforce smoothness along cortical surface
"""
batch_size, n_sources, n_time = sources.shape
# Compute pairwise differences
diff = sources.unsqueeze(2) - sources.unsqueeze(1) # [batch, n_sources, n_sources, n_time]
# Weight by inverse geodesic distance
weight = torch.exp(-self.geo_dist / torch.median(self.geo_dist)) # Closer = more similar
weight = weight.unsqueeze(0).unsqueeze(-1) # [1, n_sources, n_sources, 1]
loss = torch.mean(weight * diff**2)
return loss
def _geodesic_regularization(self, sources):
"""
Regularize based on geodesic distance along cortex
"""
# Penalize high-frequency spatial patterns
# (implementation depends on specific parcellation)
return torch.tensor(0.0) # Placeholder
Applications
- Precision Neuromedicine: Subject-specific brain activity mapping
- Real-time Neuroimaging: Online dynamics reconstruction
- Brain-Computer Interfaces: Accurate cortical state estimation
- Cognitive Neuroscience: High-resolution spatiotemporal analysis
- Clinical Diagnostics: Patient-specific brain dynamics for epilepsy, stroke
Key Features
- Individual Geometry: Subject-specific cortical surface incorporation
- End-to-end Learning: Direct sensor-to-source mapping
- Physics Constraints: Biophysically plausible reconstructions
- Geodesic Regularization: Anatomically-informed smoothness
- Multi-scale: Operates at multiple spatial resolutions
Pitfalls
- Mesh Quality: Requires high-quality individual cortical surfaces
- Computational Cost: Geodesic distance computation is expensive
- Registration: MRI-to-head coordinate alignment critical
- Head Model: BEM accuracy affects forward solution
- Training Data: Needs paired source-sensor data (simulated or intra-cranial)
Related Skills
- geometric-brain-dynamics-mapping
- eeg-visual-attention-decoding
- brain-dit-fmri-foundation-model
References
Saleh, Y., et al. (2026). A geometry aware framework enhances noninvasive
mapping of whole human brain dynamics.
arXiv preprint arXiv:2604.25592v1.