Introduction

Three decades into the human brain mapping revolution ushered by functional magnetic resonance imaging (MRI), large swathes of the neural landscape remain shrouded in darkness. Two cardinal aspects of the task are increasingly recognised to inhibit progress: the distributed, connective nature of the neural substrate1,2, and the difficulty of establishing necessity from predominantly correlative data3,4. Each aspect on its own presents formidable difficulties: characterising distributed substrates requires explicit modelling of remote interactions intractable without large-scale data and mathematical models embrittled by their complexity; establishing necessity requires disruptive evidence typically obtained naturally, through the behavioural consequences of uncontrolled focal pathological lesions confounded by their incidental—and heterogeneously distributed—characteristics5,6. Combined, these difficulties are reciprocally amplified: data of sufficient scale and quality to support complex models is especially hard to acquire in the clinical domain, and distributed patterns of pathological damage become entangled with the comparably distributed underlying patterns of neural dependence (with rare exceptions7). Yet it is precisely distributed substrates that are most in need of disruptive evidence, for the plurality of neural support makes inferences from correlative data all the harder.

Methodological innovation at the intersection of connective and disruptive mapping of brain function is therefore urgently needed, with close attention not only to the practicalities of scaling current techniques, but also to diminishing the need for data volumes that will always be hard to achieve. Here we elaborate conceptually, implement technically, and demonstrate empirically, a simple, principled approach to connective disruptive mapping of human brain function in the clinical context of direct cortical electrical stimulation (DCS).

Theoretically, the ideal approach is to register the functional consequences of transient disruption applied at single point loci, individually and in combination, across the entire brain. DCS, commonly employed as a localising clinical tool in patients undergoing evaluation for resective surgery of (typically epileptogenic) lesions, approximates this ideal arguably closer than any other available tool. Focal, transient disruption can thereby be achieved, enabling causally more robust examination of the relationship between a well-defined neural substrate and an observed, or reported, behavioural outcome8. Though clinical imperatives inevitably constrain the choice of locations and sampling density, the ability to evaluate multiple loci in each patient, dynamically, yields higher volumes of informative data than the bare number of surveyed patients suggests. The approach has already been extensively used to derive maps of functional dependence9 in surgical settings10,11, including connectivity12, but outside a formal framework that allows both focal and connective effects to be robustly quantified without dependence on predefined regions of interest.

For all its theoretical power, the use of DCS for spatial inference is complicated by its sparsity. Although it is common practice to evaluate multiple loci in each individual patient—and grid electrodes may offer locally dense coverage—comprehensive sampling at high resolution across the brain is infeasible. The traditional solution is to adopt an a priori region of interest (ROI) parcellation, and report behaviour averaged across each sampled region13,14,15. Informing regional parcellations by richer representations of neural similarity such as histology, functional, and structural connectivity (e.g. ref. 16,) increases our confidence in their fidelity, but only as far as these characteristics may reasonably be taken as indicators of functional homology: a question that can be definitively settled only by disruptive techniques itself. Moreover, this approach to anatomical inference has six defects. First, it assumes that the constituents of each ROI are homogeneous and interchangeable, a simplistically modular, “Lego” vision of the brain not sustainable on close examination17. Second, it assumes that the ROI allocation of a given stimulated locus is both certain and invariant to its distance from the ROI boundary, neither of which is plausible. Third, the resultant inference presupposes the topology it is supposed to reveal, for it is expressed in a parcellation defined before the data is even acquired. Fourth, no regional difference will register where the true functional boundary is orthogonal to the a priori one. Fifth, where a functional pattern exhibits a finer anatomical organisation than the a priori parcellation, it will be invisible through it. Sixth, both continuous and discrete spatial variations in function will appear equally abrupt.

These defects have motivated us to develop a different approach, analogous to meta-analytic mapping18,19, that enables inference to the spatial characteristics of sparsely sampled critical areas without any prior assumptions on their structure beyond a reasonable degree of local smoothness20. Consider in illustration a target anatomical domain—the dorsal medial wall, represented in 2D for simplicity—where a set of N discrete loci registered on a common grid are associated with two different deficits (Fig. 1). Though the data is in anatomical register, we cannot perform a statistical test at every point on the grid, for no point is sufficiently sampled. The conventional solution is to aggregate the observed deficits within pre-defined ROIs, and report statistics in regional terms (Fig. 1a). Framed as count regression, for example, the task is to predict the number of instances observed within an ROI (the dependent variable) given the behavioural parameters (the independent variables), iterating across ROIs. This yields a map structured by the chosen parcellation, sensitive to the correspondence between the parcellation and the underlying functional substrate. Where the two correspond poorly (Fig. 1b) the inference will be distorted or fail altogether.

Fig. 1: ROI-based vs Random Field-based Spatial Inference.
figure 1

a ROI-based spatial inference. Simulated stimulation points resulting in two hypothetical deficits—A (blue) and B (red)—are counted across pre-defined ROIs, here represented in 2D space for simplicity. A statistical test is then performed on the counts to infer the spatial distribution of function in terms determined by the ROI boundary. b Parcellation-induced mislocalisation. When there is a poor correspondence between the parcellation and the underlying functional substrate, the inference either fails completely (left) or is distorted (middle and right). c Random field-based spatial inference. Here each stimulation point is convolved with a predefined Gaussian kernel, so that each location is now supported across the entire domain, enabling the application of voxel-wise inference on a regular grid. A statistical test is then performed at each voxel to retrieve the substrates associated with the observed deficits A and B. The colourmap is the negative decimal log of each p-value. Thresholding following multiple comparisons correction is not shown here for simplicity.

Our alternative approach is to transform each sparse location from a single point to a dense spatial distribution, thereby enabling point-wise mass univariate inference on a regular grid. In the simplest, focal form, this is achieved by convolving each point with a Gaussian of predetermined width, assuming that each sampled point is drawn from an underlying spatial distribution that a random Gaussian field can approximate21. The locus of disruption is thus represented not as a point but as a continuous spatial distribution whose density gracefully captures the uncertainty of each disrupted location across the entire anatomical domain (Fig. 1c). Now framed as linear regression, the task is to predict the density observed within each voxel (the dependent variable) given the behavioural parameters (the independent variables), iterating across voxels in mass-univariate fashion. Though it may seem counterintuitive to designate location density as the dependent variable and disrupted behaviour as an independent variable, the choice here is motivated by the task of spatial inference in the context of mass-univariate analysis, and is aligned with established practice in other domains, such as voxel-based morphometry and functional imaging. For example, in voxel-based morphometry, the dependent variable is tissue concentration, and behaviour is amongst the independent variables.

In the more complex, connective form, the transformation is achieved by a probabilistic projection of the distributed connectivity of each locus, incorporating not just the uncertainty but also the network distribution of the disruption22. In both cases, this transformation allows us to use the well-established, principled approach to mass univariate spatial analysis embodied in the statistical parametric mapping (SPM) platform23.

To demonstrate the practical application of our approach, we investigate the focal and connective organisation of the medial frontal wall in the context of sensorimotor behaviours. Surveying the literature on the medial frontal cortex reveals an inferential landscape dominated by correlative studies24,25,26,27,28, disclosing a remarkable diversity of cognitive and behavioural associations, many of them conflicting. Disruptive studies are rare here for three reasons. First, stroke—the commonest source of focal lesion data—rarely involves the medial wall, and when it does, the size and morphology of the injury limit its spatial resolving power5. Second, natural lesions are never truly local: they enclose larger areas of tissue than are plausibly functionally uniform, resulting in confounding from collateral damage that is hard to remove (e.g. ref. 29). Third, non-invasive disruptive methods such as repetitive transcranial magnetic stimulation are either limited to the dorsal surface, or in reaching deeper presuppose its confounding collateral disruption along the way8. Inferences from pathological forms of focal injury are in any event complicated by plasticity and reorganisation over time30, limiting generalisability.

Here we re-examine a previously reported set of thirty-seven patients undergoing direct electrical cortical stimulation across the medial wall in the context of clinical evaluation for surgical treatment of non-lesional epilepsy20. The proximity of critical medial motor areas and the propensity for seizures to involve them justifies dense sampling of the area, achieved either by placing surface electrode grids or with multiple depth electrodes. We adapt the approach to focal voxel-wise inference pioneered by meta-analytic mapping18,19 presented in our first report20, reformulating it within SPM’s statistical framework, and extend it to disruptive connective analysis31,32 of a transient kind, introducing the notion of transient dysconnectome mapping. We speak of a ‘connectome’ because the inferred maps capture distributed anatomical relations defined by any chosen index of connectivity, and we prefix the term with ‘dys’ (rather than the usual ‘dis’) because a transient intervention typically does not disconnect a network but renders it dysfunctional. Our reformulation enhances the statistical efficiency and sensitivity of the core approach, and its extension enables us to compare focal and distributed effects, not just along the medial wall, but across remote brain regions interconnected with it, establishing a comprehensive platform for transient disruptive mapping of the human brain.

Results

Local disruptive mapping

A total of 477 disruption sites confined to the medial frontal wall were obtained from 37 patients, providing good sampling coverage across the medial wall (Fig. 2, Supplementary Table 2).

Fig. 2: Distribution of stimulation locations over the frontal medial wall.
figure 2

The mean stimulation density (smoothed by a 10 mm FWHM Gaussian kernel) is shown overlaid onto the FA template. The colourbar shows the density of electrode locations. A total of 477 stimulation points were analysed.

Positive motor responses—the most common stimulation-induced behaviour—were observed in 153 stimulations (32%) (Fig. 3, Supplementary Table 3). They were associated with disruptions in the right SMA and the precentral gyrus bilaterally. Negative motor responses were observed in 41 stimulations (8%). They were associated with disruptions in the SMA and the right middle segment of the superior frontal gyrus. Sensory responses were observed in 46 stimulations (9%). They were associated with disruptions in the left middle cingulate gyrus, and the middle segment of the precentral gyrus (bilaterally). Speech disturbances were observed in 46 stimulations (9%). They were associated with disruptions of the left superior frontal gyrus and pre-SMA. No-responses (silent disruptions) were observed in 243 stimulations (51%). They were associated with disruptions of the middle cingulate gyrus and the middle segment of the superior frontal gyrus. Results at the uncorrected threshold, are available in Supplementary Fig. 1.

Fig. 3: Local disruptive mapping of behaviour.
figure 3

For each MNI voxel, a planned t-contrast was performed. Only voxels surviving the p < 0.05 FWE-corrected threshold are shown, overlaid on the mid-sagittal plane, where higher t-statistics (brighter colour) represent a stronger association between the electrode density value and the observed behaviour. 153 stimulation points were associated with positive motor responses, 41 with negative motor responses, 46 with sensory phenomena, 46 with speech disturbances and 234 were deemed silent contacts. R right, L left, (pre-)SMA supplementary motor area.

To verify the focal results are not dependent on the predetermined smoothing kernel, we performed a sensitivity analysis using 8, 10, and 12 kernels (Supplementary Fig. 3), revealing no significant differences in the location and spread of the effects.

Dysconnectomic disruptive mapping

Positive motor responses were localized to 14 clusters (Fig. 4, Supplementary Table 4). Cortical grey matter regions included the precentral gyrus, superior parietal lobule, superior frontal gyrus, posterior insula, central operculum, and the superior occipital gyrus. Deep grey matter regions included the thalamus, ventral diencephalon, caudate nucleus, and pallidum. Cerebellar connectivity was evident in regions falling within the fronto-parietal, foot, and hand sensorimotor cerebellar networks33. At the uncorrected threshold, connectivity was also evident in the mid-brain, pons, and medulla, as can be appreciated in Supplementary Fig. 2. Negative motor responses were localized to 12 clusters, with the global maximum located in the superior parietal lobule. Cortical grey matter regions included the superior parietal lobule, SMA, superior and medial frontal gyrus, and precentral gyrus. Significant clusters were also found in the precuneus, and the cerebellar fronto-parietal and dorsal attention networks. Sensory responses were localized to 5 clusters, with the global maximum falling within the precentral gyrus. Significant areas were also obtained in the middle cingulate cortex, thalamus, the fronto-parietal cerebellar network, and the medulla. The absence of a response to disruption was associated with 7 clusters, with global maxima falling in the medial segment of the superior frontal gyrus. Other areas included the middle frontal gyrus, medial orbital gyrus, precuneus, angular gyrus and the dorsal attention cerebellar network. Note the connective maps here are normalised by distance as outlined in the methods, potentially resulting in departures from the focal maps above.

Fig. 4: Connective disruptive maps of behaviour.
figure 4

For each MNI voxel, a planned t-contrast was performed. Only voxels surviving the p < 0.05 FWE-corrected threshold are shown, overlaid on the mid-sagittal, lateral, superior, and inferior planes, where higher t-statistics (brighter colour) represent a stronger association between the connectivity value and the observed behaviour. Cerebellar subregions are labelled with reference to a priori known cortical network associations33. R right, L left, SMA supplementary motor area, SPL superior parietal lobule, SM sensorimotor.

Predictive analysis results

The foregoing maps represent population-level spatial inference to the anatomical substrate of the elicited behaviours. Their utility lies in disclosing the general anatomical organisation of the brain across the population rather than capturing individual-level variation. Population-level maps may nonetheless be useful as priors for individual localisation, in conjunction with models specifically crafted for the task. We can quantify the utility of such prior maps by examining their individual-level predictive fidelity on held-out data. Note the objective here is not to use the model as an individual-level classifier—it is not designed for that—but rather to quantify the localisation information in the resultant maps with a view to their utility as priors for individual-level modelling.

Fig. 5 shows the average ROC curve (±1 standard deviation) for 100 bootstraps of the test set for classifying each stimulation based on the spatial map inferred from the remaining data: either SPM or ROI. Predictive performance was superior with the SPM map for all behaviours. The average performance metrics in terms of accuracy, balanced accuracy, precision, true positive rate (recall), and false positive rate are provided in Supplementary Table T5. The full results for each iteration are provided in the online Supplementary Data.

Fig. 5: Average receiver operating characteristic curve (ROC) for predicting behavioural outcome from spatial maps.
figure 5

The performance of a simple classifier of the behaviour based on an SPM map (red) or an ROI map (blue), both derived from held-out data, is presented as the average of the ROC curves from 100 bootstraps of the test set. The standard deviation (± 1 SD) are plotted in shaded colours. Predictive performance was superior with the SPM map for all behaviours. Note the purpose of this analysis is not to craft an optimal behavioural classifier based on stimulation data but to quantify the comparative value of SPM vs ROI-based maps in providing spatial priors for downstream individual-level localisation.

Discussion

We have presented a framework for local and connective spatial inference with sparsely sampled focal disruptive data and applied it to transient direct cortical electrical stimulation of the medial frontal wall. Here we review the characteristics of our approach, and examine the empirical results drawn from its application to illuminating the local and distributed organisation of the medial frontal wall.

Spatial inference from sparsely sampled discrete data

Inferring a dense map from sparsely sampled disruptive data inevitably implies interpolation between unsampled regions of the target space guided by the chosen method’s inductive bias. ROI-based analysis assumes all voxels within any given region are equivalent, changing function with implausible abruptness across regions16. It yields maps that are discretised not by the data themselves, but on an assumption the inference rests on. If there is a spatial scale at which this assumption is safe, it is yet to be discovered: current parcellations of the brain are based on arbitrarily discretised continuities. By contrast, the alternative approach we pursue here merely assumes that the underlying substrate can be modelled by a random Gaussian field. The one heuristic—the width of the Gaussian with which point data is convolved—can be both empirically informed by data from other modalities and corresponds to the readily intelligible notion of chosen scale of spatial analysis.

Representing the data in Gaussian-convolved form furthermore enables the application of voxel-wise mass-univariate methods whose sensitivity and statistical efficiency makes them preferable where, as here, a given locus may licitly be evaluated independently21. SPM’s established statistical framework combines flexibility in the framing of voxel-wise statistical hypotheses with great fidelity in the context of marked noise. Where multiple loci are disrupted together, as in natural lesions, the assumptions of benignly structured local dependence are broken, and mass univariate inference is not valid, whether with dense or sparse data, and a multivariate approach becomes necessary.

Representing stimulation data as a disrupted cortical network of connections—a dysconnectome—naturally extends the evaluated neural support across the brain. Connectivity need neither be structural—functional or meta-analytic indices, for example, are equally applicable—nor unitary: multiple dysconnectome maps could theoretically be compared, within the same or parallel models. Metrics of connectivity may be transformed to focus attention on different spatial scales of interaction, such as the distance-based normalisation employed in our study to highlight long-range connectivity. Though group inference requires a common space (at least implicitly), a dysconnectome may be estimated from individual rather than template connectivity data. Indeed, where multiple data pertain to a single individual, all analysis may be conducted in the same, single native space, with the objective of making inferences about that individual alone. Dysconnectome maps are open to mass-univariate inference, with the caveat that the interactions between distributed areas are left unmodelled, just as they are when multiple regions are identified in functional imaging. A full network analysis requires an explicit model of the interactions between regions, for which a dysconnectome map may serve as an initial feature selection step. Equally, though the source data is disruptive—we are observing the alteration of cognition and behaviour by stimulation—its connective anatomical projection inevitably relies on inferred architectural commonalities.

Our approach is analogous to lesion disconnectome mapping, where a lesion is projected along the white matter tracts it encloses, yielding a probabilistic representation of disrupted white matter pathways delineated up to the boundaries of the cortical areas they connect22. The difference is that our lesions are single point loci falling entirely within grey matter, and the projection we seek is to connected grey matter areas, not to the underlying white matter. The inferred topology is then primarily in terms of grey rather than white matter anatomy, defined by any chosen mechanism of connectivity.

An obvious question is the optimal stage at which connective inference should be performed. One could estimate a dysconnectome from the regions a prior local analysis identifies as significantly associated with the behaviour at the group level. This alternative is suboptimal for two reasons. First, it precludes modelling of individual variations in the distributed substrate that individual connectivity maps could provide. More importantly, it diminishes the spatial resolution of the inference: two sites of disruption too close to be resolved within a focal map may exhibit connectivity profiles that successfully distinguish them at the dysconnectome level. One may conceive of a dysconnectome as a projection into a higher dimensional anatomical space, akin to an anatomical support vector machine, that enhances the separability of the underlying patterns of dependence.

These generalities apply to the particulars of direct cortical electrical stimulation data. Note that inter-subject variability, even after precise non-linear registration, likely dwarfs the effects of local current spread, rendering physics-informed spatial priors unnecessary34. The method of inference naturally cannot overcome the limitations of the data themselves, which here include sparse and clinically biased sampling, non-linearities of induced effects, clinically-determined variations in stimulation parameters, interactions—both enduring (e.g., coexisting structural lesions) and dynamic (e.g., after-discharges)—with the pathology that motivates intracranial study, confounding from other interventions (e.g., anti-epileptic drugs), and the relatively narrow repertoire of behaviours the clinical setting permits us to evaluate35. That the disruption is localised to a comparatively small volume, however, is a crucial advantage over natural lesions, where the distributed nature of the damage enormously complicates spatial inference, and where many of the foregoing limitations apply equally.

Clinical constraints make it inevitable that only a comparatively small number of areas can be sampled in any one patient. Just as with structural lesions, such as ischaemic stroke, where each lesion is associated with a single patient and is treated as an independent datum, so here each disruption site is treated as a single lesion and as an independent datum. Where multiple sites are drawn from the same patient, however, each observation is no longer independent, which needs to be accounted for in the statistical model. Here we employ a single-level repeated measures design, where linear variations in the mean of individual patients is absorbed by a subject-level regressor, and the non-sphericity of the error induced by within-subject correlations is modelled by whitening the design matrix36. This manoeuvre, standard in SPM, does not change the degrees of freedom of the model, unlike alternatives such as the approximate Greenhouse-Geisser correction37. A valid alternative is to adopt an approximate hierarchical, two-level design where contrasts from individual subject-specific models at the first level are evaluated by a one-sample t-test at a second level. We favour the former approach owing to its flexibility and power in the context of the imbalanced data the present clinical context tends to compel.

Special attention must be given to the heterogeneity of the neural effects of electrical stimulation. The effects at the single neuron level may be excitatory or inhibitory, varying both over time and distance, and need not correspond to effects at larger scales that inevitably depend both on individual neuronal responses and the configuration—local and remote—of physiological excitatory and inhibitory neuronal functions. Nor can a predominance of induced functional excitation or inhibition be inferred from the elicited behaviour: the interruption of movement need not be neurally inhibitory even if it is behaviourally so, and vice versa. Equally, a movement may be disrupted not because it has lost its substrate but because a competing substrate has been driven to supplant it; and may fail to manifest positively not because the substrate is not excited but because the complexity of the movement demands a pattern of excitation too elaborate for electrical stimulation to induce. These considerations are especially pertinent to the interpretation of rostral effects where negative effects or silence reign: we cannot conclude that these regions are inhibitory, only that the behaviour that depends on them may be competitive or complex. In short, electrical stimulation is correctly viewed as transient, heterogeneous disruption, of value in localisation but not necessarily the more detailed characterisation of the neural substrate.

The framework we have described is designed to reveal commonalities of functional anatomy across a population under study, not to capture heterogeneities specific to individuals. It permits principled inference to the regional organisation of the brain, from which individuals may depart in a way more flexible statistical models are required to capture. The flexibility of such—necessarily multivariate—models makes them susceptible to overfitting, especially in the low-data regimes characteristic of the domain of cortical stimulation. Our task is not to create a substitute for such models, but to provide maps that disclose the general organisation of the brain and may be useful in individual-level prediction as anatomical priors, constraining the flexibility of the predictive model. We have shown that such priors are superior when drawn from SPM rather than pre-defined ROIs weighted by exactly the same data. This observation corroborates the intuition, illustrated in Fig. 1, that data coerced into pre-defined ROIs can be expected to yield a high-fidelity map only where its boundaries naturally coincide with those of the chosen parcellation, an assumption not substantiated by the facts here, and implausible elsewhere.

Focal mapping of the medial frontal wall

Our analysis broadly replicates the maps obtained from the original permutation-based voxel-wise analysis of the same data20. A rostro-caudal organisation of behavioural complexity is observed26,38,39,40,41,42, with positive motor responses caudal to negative motor, speech, idiosyncratic phenomena (‘other’ responses), and absence of any response. In presupposing a task that is contextually interrupted, negative responses can reasonably be expected to be more complex in their condition-action association than positive ones. Speech elicited a locally distributed pattern likely reflecting the complex compositionality of the task, including relatively low-level aspects of articulation43,44.

The crudity of the behavioural labels here is a reminder anatomical inference is not to functions but to behaviours speculated to depend on them. Obviously, no function corresponds to the “silence” observed across of the vast expanse of the medial wall where no response of any kind could be elicited. The inference to be drawn is not that this region does nothing but that the behavioural tasks commonly deployed to explore its function in the clinical setting lack range and specificity.

Connectomic mapping of the medial frontal cortex with transient dysconnectomes

Projecting disrupted cortical areas to their grey-matter connections reveals a much richer picture of the underlying neural substrate. Positive motor responses are centred on the supplementary motor area, with extensive precentral gyrus involvement plausibly reflecting strong connectivity with primary motor cortex. There are extensive projections to subcortical targets including thalamus, basal ganglia, pons, medulla, and cerebellum. Negative motor responses show more rostral weighting on the medial wall, and weaker connections to precentral and deep areas. By contrast, SPL involvement was more prominent here, giving context to recent evidence for direct parietal-motor functional connectivity in the region45. Sensory responses highlight the middle cingulate gyrus in keeping with its known behavioural associations (e.g. refs. 46,47,48,49), and connectivity patterns50. The region we identify is plausibly part of the caudal cingulate premotor area involved in the multisensory orientation of the head and body in space47.

Special note should be made of cerebellar connectivity, structured by its recently described network organisation33. Here positive motor responses engage multiple areas including cerebellar nodes of the fronto-parietal, hand and foot sensorimotor, and dorsal and ventral attention networks. Negative responses are by contrast dominated by the fronto-parietal network, as befits their more complex generative context. Sensory responses overlap with a component of the cingulo-opercular network, in line with the strong association with the mid-cingulate. The association of silent responses with cerebellar connectivity within the dorsal attention network, shown to play a role in memory and attention tasks previously thought to be dominated by cortical regions51, suggests the category of tasks more likely to be eloquently modulated by disruption of the rostral medial wall. Our findings highlight the strongly structured spatial organisation of the cerebellum, inviting future research to further delineate its interaction with cortical and subcortical regions.

Conclusion

We present—and apply to the medial frontal wall—a random field-based approach to drawing spatial inferences about the focal and distributed functional organisation of the brain from sparsely sampled disruptive data. Our approach combines minimal anatomical and physiological assumptions with a principled framework for establishing disruption-behavioural associations. Applied to the medial wall, it reveals marked differences between focal and distributed maps, even in the context of relatively constrained spatial sampling, with implications for our understanding of the functional organisation of the region, and—more generally—the optimal path to integrating local and distributed information in our models of the brain.

Methods

Published data from Trevisi et al.20 were reanalysed for this study. We performed two sets of analyses. The first to derive cortical maps of focal regions critical for specific categories of sensorimotor behaviours, the second to extend these maps to connected regions across the brain within the same inferential framework.

Participants

Retrospective data from 147 consecutive drug-resistant focal epilepsy patients undergoing intracranial recording between January 2008 and June 2015 at the National Hospital for Neurology and Neurosurgery, London as part of clinical assessment prior to epilepsy surgery were screened. The study was approved by the hospital as a retrospective evaluation of routine clinical practice.

Thirty-seven patients (28 males, 9 females, aged 19–68 years, mean = 33.86 years, SD = 11.11 years) were identified to have at least one stimulation in the medial frontal region, spatially confined to the medial wall area dorsal to the corpus callosum and rostral to the caudal bank of the marginal sulcus. Five patients had minor lesions near the supplementary motor area (SMA), and one had a lesioned SMA/paracentral lobule as evident on MRI. Table 1 summarises the patients’ demographic and clinical characteristics [based on Trevisi et al.20]. Further details are provided in Supplementary Table 1.

Table 1 Demographic and clinical characteristics.

Direct cortical stimulation procedures

As previously described20, depth electrodes were implanted in 19 patients using a frameless stereoelectroencephalography (SEEG) technique52. In the remaining 18 patients, craniotomy was performed for the placement of strips and/or grids with or without freehand insertion of additional depth electrodes. In 14 patients, intracranial recording was performed in the right hemisphere, in 18 patients in the left hemisphere, and in 5 patients bilaterally. Hemispheric dominance for language was inferred from fMRI data, not re-examined here. In 19 (51%) cases, the recording was in the dominant hemisphere, whereas in 18 (49%) cases, the electrodes were in the nondominant hemisphere. Bilateral language dominance was noted in five patients. The location of the electrodes was confirmed for all patients by post-implantation CT studies.

A clinical epileptologist and a physiologist performed one or more sessions of DCS during simultaneous video-EEG recording. Stimulations were typically performed after ictal recordings when patients were back on their baseline antiepileptic medication. Bipolar or monopolar stimulation trains were delivered with biphasic rectangular pulses of AC-current at 50 Hz, with a pulse width of 500 µs and a maximum duration of 5 s. The intensity was gradually increased from 0.5 to 7 mA in increments of 0.5–1 mA until the occurrence of a clinical sign or until after-discharges were detected on EEG monitoring20,53. Full details on the electrical stimulation intensities for each of the behavioural categories are provided in Trevisi et al.20. Stimulations accompanied a stereotyped set of test actions—rest, Barré and/or Mingazzini test, repeated movements of the upper and lower limbs, and during counting, reading, or repetitive monosyllabic verbalisation—as described in detail elsewhere20.

Behavioural analysis

Three clinicians classified the observed behavioural responses as ‘positive motor’, ‘negative motor’ or ‘speech disturbances’. Positive motor responses included involuntary, typically tonic or clonic, movements of the eye, head, limb or trunk. Negative motor responses included slowed or inhibited movement relative to experimentally specified movements, such as the inability to maintain prescribed postures. Speech disturbances included speech arrest, alteration in rhythm, involuntary speech, and hesitation. Live and post-hoc classifications (using video and audio telemetry recordings) were made by stimulating and attending clinicians. A contact was deemed silent if no response was obtained at the maximum stimulation intensity of 7 mA (‘no response’). Responses after a seizure or after-discharges were excluded from the analysis.

Patients were also asked to report any evoked somatosensory responses such as cutaneous paraesthesias (tingling, touch, heat, and pain). Responses that were neither sensory nor motor, such as a reported urge to move or speak, or reported perception of motion without observed movement, were elicited only three times across the entire dataset and were therefore not modelled. Note that since each locus was evaluated with multiple tasks, more than one class of response may be associated with it: the classification is not anatomically exclusive.

Imaging data acquisition and processing

Preoperative structural T1-weighted imaging with an isotropic resolution of ~1 mm was acquired on a 3 T magnetic resonance imaging (MRI) scanner. After implantation, non-contrasted structural whole-head CT scans with a resolution of 0.43 × 0.43 × 1.2 mm (SOMATA Definition 128-slice, Siemens Healthcare GmbH, Erlangen, Germany) were obtained to confirm the location of electrode contacts. All image processing was performed using SPM12 (http://www.fil.ion.ucl.ac.uk/spm/).

To facilitate group analysis, electrode locations were manually extracted from the CT and non-linearly transformed into the Montreal Neurological Institute (MNI) space template as described in Jha et al.7. In brief, for each patient, a rigid body co-registration to the standard SPM12 tissue probability map was performed for both preoperative T1-weighted MRI and postoperative CT, based on normalised mutual information with adjustment by Procrustes analysis, weighted by white and grey matter compartments. The algorithm brought each scan into an approximately rigid register with the MNI template, making subsequent transformations more robust. The standard co-registration algorithm in SPM12 was then applied to co-register each CT scan with its coupled T1-weighted MRI, which enabled automatic replication of each subsequent transformation of the T1-weighted images with their corresponding CT scans. Standard segmentation and normalisation routines with default parameters were then applied to the T1-weighted images to create segmented images in native space for each of the six standard tissue classes, combined with a set of non-linear parameters, to transform the resulting segments into MNI space. These parameters were then used to transform the white matter and grey matter compartments of each T1-weighted image and the corresponding CT scan into normalised MNI space. The location of the electrodes in MNI space was determined by displaying the normalised T1-weighted and CT images together in triplanar view using the SPM12’s ‘check registration’ module. The location of the centre of each electrode was visually judged by two independent observers to lie within the grey matter of the medial wall. By locating the electrode contacts in MNI space after normalisation, the potential bias due to anatomical differences between subjects was minimised during manual labelling.

Local disruptive mapping

A total of 477 stimulation locations were extracted, covering a larger region than the original study. For each location a corresponding image (1.5 × 1.5 × 1.5 mm sampled) was generated with intensity of zero at all locations except the stimulated location where the intensity was one. Each image was convolved with a 3D Gaussian 10 mm full-width-half-maximum (FWHM) kernel (truncated at 90% mass) to enable modelling of spatial uncertainty in the location of the stimulation and approximate the local distribution of focally induced disruption. The resultant image contained a single Gaussian located at the stimulation point, representing location uncertainty as the density of this spatial distribution across all voxels. This approach facilitated group analysis of sparse data accounting for between-subject variation in functional-anatomical relationships not captured by anatomical registrations, analogous to the approach used in meta-analytic modelling of functional activation data18,19. The kernel size used here was guided by empirical studies on spatial uncertainty modelling of functional neuroanatomical data12. Manipulating the kernel width from 4 to 16 mm in increments of 2 mm, Trevisi et al.20 reported similar results, suggesting that the choice of the kernel size is not a critical step in the analysis. We do not claim that our choice of kernel width is optimal or generally prescribable, though it can be—as here—empirically informed by data from other modalities. Note also that the kernel size is dominated by plausible inter-subject variability rather than the comparatively much smaller scale of current spread34.

Trevisi et al. focus on the rostro-caudal organisation of the medial wall motivated them to collapse the data across other planes20. Here the enhanced efficiency of our approach allowed us to investigate bilateral effects. Data were masked by applying a threshold where electrode density was >0.00001 to exclude areas with poor sampling. The subsequent mask was confined to the frontal medial wall, and extended laterally 22 mm to encompass its depths. For each behavioural condition of interest, stimulation images were entered into a voxel-wise repeated-measures general linear model with electrode density as the dependent variable and subject and the binary behavioural effect as the independent variables. Within-subject non-sphericity of errors was accounted for using standard procedures36. A planned one-tailed voxel-wise t-test of each behavioural condition was performed and thresholded at p < 0.05 FWE (peak voxel) to account for multiple comparisons.

Connective disruptive mapping

We used large-scale high-resolution diffusion tensor imaging from the Human Connectome Project (HCP) to derive white matter connectivity matrices54,55. The imaging acquisition protocols for the HCP are described elsewhere54. In total, data from 945 participants were deemed suitable for analysis. We used FSL pre-processed diffusion data supplied by the HCP working group54. In brief, this pre-processing included b0 signal intensity normalization across the six-diffusion series, and correction for echoplanar imaging distortion, eddy current and subject motion distortion, and gradient nonlinearities. Registration of the diffusion images to the native T1-weighted structural space enabled non-linear registration to 2 × 2 × 2 mm sampled isotropic MNI template space via a deformation field derived with FNIRT from each individual’s T1-weighted image. This included BEDPOSTX processing with the default deconvolution model using sticks with a range of diffusivities.

Probabilistic tractography was applied to the diffusion data to derive local fibre orientation information. Tractography was performed using the GPU Bayesian implementation of probtrackx256, to derive a network representation of white matter structural connections across all grey matter voxels at 2 × 2 × 2 mm resolution for each of the 945 participants. An image resolution lower than the source acquistion was chosen owing to the considerable computational challenge of estimating an network of this size. Our processing parameters included passing a grey matter mask, 5000 samples, a curvature threshold of 0.2, 2000 steps with a steplength of 0.5 mm, and a subsidiary fibre volume fraction threshold of 0.01, with normalisation by the participant waytotal. The waytotal is the total number of generated tracts that satisfy the inclusion/exclusion mask criteria: normalising by it scales the estimated values to the local tract density, enabling better accounting of connectivity variations with distance. The probabilistic tractography approach and implementation are described elsewhere56,57.

Having derived a white matter tractography network of each grey matter voxel, we averaged the streamline samples across all 945 patients to yield a large adjacency matrix, which could be incorporated into an undirected weighted graph. The graph comprised 125760 individual grey matter voxels, with 7907725920 unique edges weighted by the mean normalised-streamline value, subsequently used to infer the strength of the structural connection between grey matter voxels.

Note that optimal approach to deriving white matter, amongst other, connectivity maps is a subject of intense study58,59: for our purposes a widely used exemplar is sufficient. Any alternative may be substituted, including subject-specific maps derived from individual tractographic or resting state imaging. Note also that the waytotal distance normalisation employed here deliberately magnifies the remote effects it is the objective of this approach to reveal. Others may choose to forego this step, or to explore its effects on the downstream statistics.

Given the set of stimulation coordinates, we used this structural connectome to reconstruct brain maps depicting connection strength from a stimulation seed point to all other grey matter voxels. The connection strength (edge weight) from the stimulated voxel to each other grey matter voxel was rendered as an image volume—one for each stimulation location—for subsequent analysis. We used intensity clamping outside 0.1 and 99.9% of the intensity cumulative density to eliminate the influence of presumably spurious extreme values.

As the resulting maps are already dense (in contrast with the sparse focal disruption maps), for the dysconnectomic analysis we applied a smaller smoothing kernel of 6 mm FWHM prior to assessing association for each behavioural outcome with the same statistical design as used above for local disruption mapping. Again, a voxel-wise repeated-measures general linear model with subject and the binary behavioural effect as factors was used, and non-sphericity of errors was accounted for36. Planned one-tailed voxel-wise t-tests were performed and thresholded at p < 0.05 FWE (peak voxel)—now revealing dysconnectome maps.

Visualisation

Visualisation was done using the SurfIce toolbox (https://www.nitrc.org/projects/surfice). We used FSL’s HCP1065 standard-space FA atlas to generate a background FA template60 converted into a mesh using SurfIce’s volume-to-mesh function.

The local disruption and normalised disconnectome maps were visualised at the FWE corrected threshold in the mid-sagittal, lateral, superior, and inferior views. The mid-sagittal laterality was determined by the x coordinated value of global maxima of each map, so that for x < 0 the left hemisphere is shown, whereas for x ≥ 0 the right hemisphere is shown. The nearest grey matter location to the maxima of each statistical cluster was determined using SPM12’s neuromorphometrics grey matter atlas for cortical regions. For cerebellar regions, we used a well-established functional parcellation into distinct networks based on resting state connectivity33. To illustrate the extent of the effects, the local disruption and normalised dysconnectome maps were also visualised at a lower statistical threshold (p = 0.001 uncorrected, t = 3.11) and are shown in Supplementary Figs. S1 and S2 respectively.

Predictive analysis

Ours is a method for population-level spatial inference, not individual-level prediction. Inferred spatial maps may nonetheless be used as spatial priors for individual predictive models, where population anatomy may helpfully constrain the space of individual variability. Here we therefore evaluate the predictive performance of SPM-derived maps for each of the four classes of elicited behaviour, with conventional ROI-derived maps for comparison.

We obtained bootstrap estimates of predictive performance for each behaviour, resampling the data across 100 mapping and testing sets, with an 80–20 split. Each mapping subset was used derive a spatial map of significant associations with either SPM as described above, or Fisher’s exact test applied to an a priori ROI parcellation. At test time, the corresponding map—SPM or ROI—was used to assign a “weight” to each stimulation locus in the test set, indicating the strength of its association with the corresponding behaviour as captured by the statistical test on which the map is based. For the SPM maps, this weight was derived from the statistic of each test voxel; for the ROI maps it was derived from the statistic of the ROI enclosing each test voxel. In every case, the voxel-level weight was used as a single predictor variable for classifying each behaviour, independently across separate models. Performance was quantified in terms of standard classification metrics—accuracy, balanced accuracy, precision, true positive rate (recall), and false positive rate—with estimates of uncertainty from the boostrap.

Specifically, for the ROI models, we used Glasser’s grey matter atlas16 to create a priori parcellation of the region of the medial-frontal wall sampled in our study. The binarized sum image of all smoothed stimulation points (thresholded at >0.0001) was multiplied voxel-wise with the original whole-brain atlas to identify 81 unique sub-regions intersecting with the sampled region. Each stimulation locus was assigned to the ROI that enclosed it, and Fisher’s exact test was run for each behaviour and each ROI on the data from each training-test split, yielding a set of regional statistics (odds ratios and FDR-corrected asymptotic p-value) quantifying the association between them. For each behaviour separately, test loci were individually labelled as predicting the behaviour or not dependent on whether they fell within ROIs significantly associated with that behaviour. The resultant confusion matrix was used to derive performance metrics as above.

For the SPM models, we replicated the approach described in the preceding sections with only subsets of the data, employing the same bootstrap splits as for the ROI models. The resultant continuous maps were thresholded at the corresponding FWE critical t-statistic, yielding binary maps that were then used as predictors for the held-out test set of each fold. As for the ROI models, for each behaviour separately, test loci were individually labelled as predicting the behaviour or not dependent on whether they fell within the region significantly associated with that behaviour. The resultant confusion matrix was used to derive performance metrics as above.

Finally, we generated average ROC curves and AUC scores with unthresholded predictors—t-values for SPM maps and ORs values for ROI maps—as input, capturing the average performance of the models across folds and their standard deviation.

Statistics and reproducibility

A total of 477 unique stimulation points were analysed using SPM12 (see full list in Supplementary Table 2), following standard procedures. The analysis was repeated with random behavioural vectors to verify the results can be attributed to the independent variable. Predictive analysis was performed in Python, where standard classification metrices were computed across all 100 bootstrap samples (raw data is available in the online version of this paper). Code is available from the first author upon reasonable request.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.