Stimulation mapping and whole-brain modeling reveal gradients of excitability and recurrence in cortical networks

Abstract

The human brain exhibits a modular and hierarchical structure, spanning low-order sensorimotor to high-order cognitive/affective systems. What is the mechanistic significance of this organization for brain dynamics and information processing properties? We investigated this question using rare simultaneous multimodal electrophysiology (stereotactic and scalp electroencephalography - EEG) recordings in 36 patients with drug-resistant focal epilepsy during presurgical intracerebral electrical stimulation (iES) (323 stimulation sessions). Our analyses revealed an anatomical gradient of excitability across the cortex, with stronger iES-evoked EEG responses in high-order compared to low-order regions. Mathematical modeling further showed that this variation in excitability levels results from a differential dependence on recurrent feedback from non-stimulated regions across the anatomical hierarchy, and could be extinguished by suppressing those connections in-silico. High-order brain regions/networks thus show an activity pattern characterized by more inter-network functional integration than low-order ones, which manifests as a spatial gradient of excitability that is emergent from, and causally dependent on, the underlying hierarchical network structure. These findings offer new insights into how hierarchical brain organization influences cognitive functions and could inform strategies for targeted neuromodulation therapies.

Fig. 1. Studying Resting-State Network (RSN) input processing strategies and the role of recurrent feedback with computational brain network models.

Shown here is a schematic of the hypotheses, methodology, and general conceptual framework of the present work. A Intracerebral electrical stimulation (iES) applied to an intracortical target region generates an early (~20-30 ms) response (evoked-related potential (ERP) waveform component) at high-density scalp electroencephalography (hd-EEG) channels sensitive to that region and its immediate neighbors (red arrows). This also appears in more distal connected regions after a short delay due to axonal conduction and polysynaptic transmission. Subsequent second (~60–80 ms) and third (~140–200 ms) late evoked components are frequently observed (blue arrows). After identifying the stimulated network in this way, we aim to determine the extent to which this second component relies on intrinsic network activity versus recurrent whole-brain feedback. B Schematic of the hierarchical spatial layout of canonical RSNs as demonstrated in Margulies and colleagues, spanning low-order networks showing greater functional segregation to high-order networks showing greater functional integration. Networks are distributed based on their position along the first principal gradient. The stimulation sites are distributed across different levels of this gradient. C Schematic of virtual dissection methodology and key hypotheses tested. We first fit personalized connectome-based computational models of iES-evoked responses to the hd-EEG time series, for each patient and stimulation location. Then, precisely timed communication interruptions (virtual dissections) were introduced to the fitted models, and the resulting changes in the iES-evoked propagation pattern were evaluated. We hypothesized that lesioning would lead to activity suppression (C, right side) in high-order but not low-order networks.

Fig. 2. Empirical high-density electroencephalography (hd-EEG) and stereotactic electroencephalography (sEEG) signals show larger global activation patterns for high-order than low-order brain networks.

A The histogram illustrates the distance in centimeters between the electrode’s centroid delivering the electrical stimulus and the center of the nearest Schaefer’s parcel. The results indicate a high level of spatial precision, with 97.2% of sessions showing distances of less than 1 cm. B Global mean field power (GMFP) of hd-EEG averaged across all 36 subjects and 323 sessions, revealing three consistent response peaks/clusters within strict confidence intervals at ~40 ms, ~80 ms, and ~370 ms, consistent with prior electrophysiological research. C GMFP of every stimulated Resting-State Network (RSN) for hd-EEG (top row) and sEEG (bottom row). The bar plot of the normalized area under the curve (AUC) of the three clusters revealed a significantly stronger global activation pattern when the stimulus targeted high-order networks, such as the Default mode network (DMN) and Frontoparietal Network (FPN), particularly for the late evoked responses (third cluster at ~370 ms). Data are presented as mean values ± standard error of the mean (SEM) (error bars), with individual subject data points overlaid (36 independent subjects, 323 stimulation sessions). In the GMFP time course plots, shaded areas represent ±SEM around the mean. Notably, this trend aligns with the “principal gradient” hierarchy reported in the functional magnetic resonance imaging (fMRI) literature, which describes a general pattern from low-order to high-order regions.

Fig. 3. Removing recurrent connections to isolate the stimulated network suppresses late evoked potentials for high-order networks.

A Global mean field power (GMFP) for every stimulated network for model-generated high-density electroencephalography (hd-EEG) data run with both the intact (continuous line) and disconnected (dashed line) structural connectome. Findings show a more pronounced decrease in evoked late responses for high-order networks (LN Limbic Network, SN Salience Network, DAN Dorsal attention network, FPN Frontoparietal Network, DMN Default mode network). B Area under the curve (AUC) differences comparing the simulation run with the intact versus the lesioned structural connectome. The bar plot shows differences across three time windows (1st response: 0−37 ms, 2nd response: 37–78 ms, 3rd response: 78–373 ms). Data are presented as mean values ± standard error of the mean (SEM) (error bars), with individual subject data points overlaid (36 independent subjects, 323 stimulation sessions). A significant reduction in the AUC was found for late responses (78−373 ms) of high-order networks (LN, SN, DAN, FPN, and DMN) compared to low-order networks (Visual Network [VN] and Somatomotor Network [SMN]), indicated by asterisks (*P < 0.05). C Demonstration of the network recurrence-based theory for two representative sessions. Simulations of evoked dynamics are run using the intact (left) and lesioned (right) anatomical connectome. In the latter case, the connections were removed to isolate the stimulated networks for SMN (top) and DMN (bottom). In the case of the low-order network, this virtual dissection does not significantly impact the evoked potentials, while for the high-order network, a substantial reduction or disappearance of evoked components was observed.

Fig. 4. Methodological workflow for characterizing the stimulated network and performing subject-specific connectome-based neurophysiological modeling of evoked potentials.

A Simultaneous stereotactic electroencephalography (sEEG) and scalp high-density electroencephalography (hd-EEG) signals were recorded. The black triangle and dashed vertical line indicate the time at which intracerebral electrical stimulation (iES) was delivered. For further details on the methodology and data preprocessing please refer to refs. ^,. B To pinpoint the brain network where the stimulus was delivered, we employed the Schaefer atlas, which divides the brain into 1000 regions across seven distinct Resting-State Networks (RSNs): Visual Network, Somatomotor Network, Limbic Network, Dorsal attention network, Ventral Attention Network, Frontoparietal Network and Default Mode Network. Subsequently, we identified the parcellation region that overlapped with the intracerebral electrode responsible for delivering the stimulus, ultimately enabling us to determine the stimulated network. C To model individual stimulus-evoked time series, the Jansen-Rit model, a neural mass model comprising pyramidal, excitatory interneuron, and inhibitory interneuron populations, was embedded in every parcel of the lower-resolution 200-region Schaefer atlas for simulating and fitting neural activity time series. The connectivity between regions was modeled using diffusion-weighted magnetic resonance imaging (MRI) tractography computed from a sample of healthy young individuals from the Human Connectome Project (HCP) Dataset, and then averaged to give a grand-mean anatomical connectome. The iES-induced depolarization of the resting membrane potential was modeled by a perturbing voltage offset to the mean membrane potential of the excitatory interneuron population. Next, a lead field matrix was employed to project the time series from the cortical surface parcels into EEG channel space, resulting in the generation of simulated scalp hd-EEG measurements. The quality of fit (loss) was quantified by calculating the cosine similarity between the simulated and empirical stimulus-evoked time series. Optimization of model parameters was accomplished by leveraging the autodiff-computed gradient between the objective function and the model parameters, employing the ADAM algorithm. Ultimately, the optimized model parameters were utilized to generate the fitted, simulated (optimized) stimulus-evoked hd-EEG activity.