Imaging the effective networks associated with cortical function through intracranial high-frequency stimulation.

Direct electrical stimulation (DES) is considered to be the gold standard for mapping cortical function. A careful mapping of the eloquent cortex is key to successful resective or ablative surgeries, with a minimal postoperative deficit, for treatment of drug-resistant epilepsy. There is accumulating evidence suggesting that not only local, but also remote activations play an equally important role in evoking clinical effects. By introducing a new intracranial stimulation paradigm and signal analysis methodology allowing to disambiguate EEG responses from stimulation artifacts we highlight the spatial extent of the networks associated with clinical effects. Our study includes 26 patients that underwent stereoelectroencephalographic investigations for drug-resistant epilepsy, having 337 depth electrodes with 4,351 contacts sampling most brain structures. The routine high-frequency electrical stimulation protocol for eloquent cortex mapping was altered in a subtle way, by alternating the polarity of the biphasic pulses in a train, causing the splitting the spectral lines of the artifactual components, exposing the underlying tissue response. By performing a frequency-domain analysis of the EEG responses during DES we were able to capture remote activations and highlight the effect's network. By using standard intersubject averaging and a fine granularity HCP-MMP parcellation, we were able to create local and distant connectivity maps for 614 stimulations evoking specific clinical effects. The clinical value of such maps is not only for a better understanding of the extent of the effects' networks guiding the invasive exploration, but also for understanding the spatial patterns of seizure propagation given the timeline of the seizure semiology.

direct electrical stimulation

functional magnetic resonance imaging

high‐frequency stimulation

stereo‐electroencephalography

INTRODUCTION

A significant amount of research has been dedicated to studying resting‐state networks in the human brain using a variety of modalities, including functional magnetic resonance imaging (fMRI), electroencephalography (EEG), and magnetoencephalography (MEG). Task‐related functional connectivity has been extensively studied using fMRI, however this method has a low time resolution. By contrast, EEG/MEG has an excellent time resolution, at the expense of a lower spatial resolution. Task‐related connectivity has been so far studied using noncausal connectivity estimates. A method providing unambiguous causal relationships between brain areas refers to actively perturbing one area of the brain using a form of stimulation (direct electrical stimulation—DES, transcranial magnetic stimulation—TMS) and measuring the effects of this perturbation on EEG signals. This approach has been used successfully to map resting‐state brain connectivity (Matsumoto et al., 2004) by measuring cortico‐cortical evoked potentials (CCEP) evoked by single‐pulse electrical stimulation (SPES; Valentín et al., 2005). Recently, a spectral‐domain analysis of the analysis shows that responses to SPES can also be used to determine interlobar effective connectivity supporting language‐related brain function (Sonoda et al., 2021). A traditional method for mapping brain function is the intracranial high‐frequency electrical stimulation (Borchers, Himmelbach, Logothetis, & Karnath, 2012; Selimbeyoglu & Parvizi, 2010). Stereoelectroencephalographic (SEEG) presurgical investigations in patients with drug‐resistant epilepsy offer a fantastic window of opportunity allowing to study cortical function, as high‐frequency stimulation (HFS) is the routine protocol for mapping eloquent cortex (Munari & Bancaud, 1987). One of the misconceptions related to DES is that it allows drawing unequivocal conclusions about the role of stimulated brain areas (Borchers et al., 2012). It is rather considered that the stimulated point is just an input gate into a large distributed network (David, Bastin, Chabardès, Minotti, & Kahane, 2010). Evidencing the networks activated during high‐frequency DES evoking symptoms is subject to several methodological challenges, the most important one being perhaps the huge stimulation artifacts obliterating the recorded EEG signals, regardless of the distance from the stimulation site. A widely used artifact reduction method that consists in alternating the polarity (AP) of the pulses in a stimulation train has been successfully used in the past for resolving evoked compound action potentials (eCAPs) in the cochlear nerve (Hughes, Goehring, & Baudhuin, 2017) and more recently in intraoperative peripheral nerve monitoring (Wu et al., 2021). The canonical AP signal analysis is performed in time domain, by averaging the responses to positive and negative polarity pulses. We propose an original analysis of the responses in the frequency‐domain that takes into account the difference between nonlinear response of the AP stimulated neural tissue and the linearity of the stimulation artifact resulting in nonoverlapping spectral lines to fully separate the evoked responses from the stimulation artifact, without blanking any part of the signal. Given that cortical responses to HFS trains are essentially steady‐state evoked potentials, this is in line with the frequency‐domain analysis normally performed in steady‐state visually evoked potentials (SSVEP) studies (Norcia, Appelbaum, Ales, Cottereau, & Rossion, 2015). All the concepts we will be using related to spectral lines, their harmonics, linear versus nonlinear responses, modulation, phase/polarity reversal have a direct correspondence in the SSVEP conceptual and methodological framework. Using such an approach allowing to disambiguate linear artifactual components from nonlinear steady‐state electrically evoked potentials (SSEEP), we show that it is possible to reveal the brain networks associated with symptoms elicited by DES.

The purpose of epilepsy surgery is to provide seizure freedom without adverse events or permanent neurological deficits. To fulfill this goal, electrical cortical stimulation has played a pivotal role in mapping cortical functions and determining indispensable eloquent areas (Kahane, Minotti, Hoffmann, Lachaux, & Ryvlin, 2003). However, the cortical network organization generating a particular clinical effect during stimulation is yet to be determined, as no direct or effective measure is currently available to quantify the stimulation‐evoked response. Several reports characterize cortical interdependencies during clinical effects elicited by stimulation (Bartolomei et al., 2012; Koubeissi, Bartolomei, Beltagy, & Picard, 2014) but they are using a functional estimation of this interaction, not an effective connectivity approach. Mapping effects' networks would help us better understand and foresee postsurgical neurological deficit and would also enable us to plan accurate SEEG strategies by linking seizure semiology with brain's functional networks.

In this study, we aim to test the hypothesis that the symptoms evoked by intracranial HFS are not only the result of local activation of the cortex (canonical view), but also the result of the activation of a broader functional network.

PATIENTS AND METHODS

Patient selection

We selected 26 patients explored by stereo‐electroencephalography (SEEG) in the Emergency University Hospital Bucharest between 2017 and 2021 (Table 1). Of the 35 consecutive patients where the study's methodology was applied, three were excluded as they had modified anatomy or prior resections, precluding accurate pial surface reconstruction and template registration, three were excluded due to a significant number of noisy or artifacted contacts and one was excluded as it was a re‐implantation with a small number of electrodes (5). No additional exclusion criteria based on gender, ethnicity were applied. All patients were diagnosed with focal drug‐resistant epilepsy of structural etiology and were admitted as possible surgical candidates. They initially underwent phase one noninvasive presurgical evaluation that included patient and family history followed by video‐electroencephalography and neurocognitive evaluation. Consequently, each patient underwent 1.5 or 3 T magnetic resonance imaging (MRI) (isotropic 3D T1, axial and coronal FLAIR and T2, T2*or SWI), and functional imaging (interictal FDG‐PET‐CT scan). For all patients, invasive exploration using SEEG was considered necessary to delineate the epileptogenic zone (Kahane & Landre, 2008; Munari & Bancaud, 1987), map the functional cortex and to determine the limits of the resection (Cardinale et al., 2013; Isnard et al., 2018; Jayakar et al., 2016; Kahane et al., 2003; Munari et al., 1994).

TABLE 1

List of patients included in this study, gender, age at implantation, age at epilepsy onset, lateralization, implantation hypothesis, number of implanted electrodes, and the total number of intracranial contacts

Patient	Sex	Age	Epilepsy onset	Lateralization	Epilepsy	Electrodes	Contacts
1	M	39	16	L	Frontal	16	209
2	M	22	10	B	Temporal	19	232
3	M	39	13	R	Parietal	11	171
4	M	31	26	R	Frontal	8	98
5	M	41	3	L	Insular	17	209
6	M	38	34	L	Opercular	9	92
7	F	19	2	L	Parietal	11	154
8	M	20	2	L	Frontal	14	174
9	F	24	9	L	Insular	14	163
10	M	28	26	R	Parietal	17	190
11	F	19	2	L	Parietal	13	155
12	F	17	13	L	Insular	12	168
13	F	26	17	R	Temporal‐occipital	18	258
14	M	32	13	R	Temporal‐occipital	15	214
15	M	19	14	R	Temporal‐insular	15	213
16	F	26	17	R	Frontal	9	100
17	F	22	3	L	Parietal	15	191
18	M	47	46	L	Temporal	10	124
19	M	38	27	R	Insular	14	172
20	M	27	4	L	Insular	10	145
21	M	26	23	R	Temporal	11	152
22	F	26	5	L	Rolandic operculum	10	135
23	F	39	29	R	Temporal	9	129
24	M	24	7	L	Insular	13	127
25	M	18	10	R	Frontal	13	189
26	M	31	1	B	Temporal	14	187

Data collection and analysis workflow

The workflow of the analysis is summarized in Figure 1, individual steps being described in detail in the following paragraphs. In summary, responses during intracranial high‐frequency electrical stimulation were analyzed and effective connectivity was calculated between the stimulated site and response sites. The locations of the recording sites, grouped by the clinical effects they evoked, were pooled at the level of the entire patient group using standard intersubject registration (Fischl, Sereno, Tootell, & Dale, 1999) and assigned a label based on the HCP‐MMP parcellation (Glasser et al., 2016). Average responses of sites having the same label were used to create per‐symptom color‐coded activation maps.

FIGURE 1

Analysis workflow that describes the process of creating per‐symptom activation maps from HFS stimulation‐evoked EEG responses and clinical effects

Invasive exploration using SEEG and data acquisition

SEEG exploration was performed using depth electrodes (Dixi, Chaudefontaine, France and AdTech, Oak Creek, WI) with 8–18 contacts per electrode, 2 mm contact length, 3.5 mm contact spacing, and 0.8 mm diameter. Multiple electrodes were placed following an individual hypothesis. Electrodes were placed intracranially using the Leksell stereotactic frame (Elekta AB, Stockholm, Sweden or the microTargeting™ Multi‐Oblique Epilepsy STarFix Platform (FHC Inc, Bowdoin, ME) (Balanescu et al., 2014; D'Agostino, Kanter, Song, & Aronson, 2019; Dewan et al., 2018; Yu, Pistol, Franklin, & Barborica, 2018). To determine the exact location of each electrode and contact, the postimplantation CT scan was co‐registered with the preimplantation MRI. Video‐SEEG recordings were performed in chronic conditions for 7–14 days at the University Emergency Hospital Bucharest. Up to 128 out of a maximum of 258 implanted contacts were continuously recorded at a sample rate of 4,096 Hz using a 64‐channel Nicolet Wireless Amplifier or two 128‐channel XLTek Quantum Amplifiers (Natus Neuro, Middleton, WI) bridged together.

Electrical stimulation

Multiple brain DES protocols were carried out as part of the standard presurgical assessment in our center: low‐frequency 1 Hz (Munari et al., 1993), high‐frequency (HFS) 40‐50 Hz (Bernier et al., 1990; Trebuchon & Chauvel, 2016) as well as SPES (Donos, Mălîia, et al., 2016; Donos, Mîndruţă, et al., 2016; Valentín et al., 2002). Constant‐current (up to 3 mA for all protocols except SPES, where the maximum current was 5 mA) biphasic square pulses having a duration of 3 ms for low‐frequency and SPES stimulation and 1 ms for HFS were applied through pairs of adjacent contacts using a programmable clinical stimulator capable of generating arbitrary waveforms (Guideline4000LP+, FHC, Bowdoin, ME). The standard duration of the HFS train was 5 s, but some stimulations were terminated early if major clinical effects or ictal symptoms were elicited. Those prematurely terminated trials were excluded from our analysis. Specific to this study is the application of a modified HFS protocol where the polarity of the biphasic pulses is being alternated (Figure 2a). This enables the separation of the tissue response from the stimulation artifact, as it will be described in detail in Section 2.6.

FIGURE 2

(a) The regular biphasic stimulation waveform , the modulating function and the alternating polarity waveform ; (b) the Fourier spectrum of the waveforms and , where and are the absolute values of the Fourier coefficients; (c) the Fourier spectrum of the modulated waveform that has null components (Fourier coefficients ) at the frequency of the original biphasic train and its harmonics

We have analyzed the EEG responses during the application of high‐frequency electrical stimulation with a frequency of 50 Hz (patients 1–12) or 43.2 Hz (patients 13–26) at current levels large enough to evoke a clinical effect. Stimulation frequency was decreased to 43.2 Hz starting patient 13 to reduce the confounding effect of the response at 50 Hz and the interference from 50 Hz power lines in Europe. The particular 43.2 Hz frequency was imposed by the hardware characteristics of the stimulator (using an integer number of 1,111 samples per cycle at a sampling frequency of 48,000 Hz), the ability to perform a bandpass filtering of the response on the stimulation frequency with a pass band of at least 4 Hz, while significantly attenuating the neighboring 50 Hz line frequency, and avoidance of the overlap not only between the two fundamental frequencies, but also of harmonics.

Clinical effects

Multiple stimulation trials were performed, gradually increasing the current intensity from 0.1 to 3 mA, usually 0.25 mA steps. Wider steps of 0.5 or even 1 mA may have been used if repeated anterior stimulations did not elicit any effect. In addition, narrow steps of 0.1 mA were used when a clinical effect was obtained at low intensities or when there was a need to choose between simultaneous multimodal clinical effects. The upper current limit was set by the presence of clinical or electrical response (after discharges). If none of these were obtained, the upper limit was set by default to 3 mA, which generates a maximum charge‐per‐density allowed in our set‐up. Patients were instructed to report any psychological or physical changes they experienced during or after each stimulation. Furthermore, patients were also engaged in different tasks during the stimulation session (e.g., repetitive movements of the limbs to test for motor deficit, reading, counting to test for language or attentional deficits, etc.). Tasks were selected based on the region stimulated. For example, for stimulations included in this study, if the contacts were placed in the proximity of the motor area (primary or secondary), patients were asked to perform a repetitive movement and to count. Tasks have no interpatient variability. The electrical stimulation was delivered for 1–3 s after the clinical testing started. In many situations, the stimulation was repeated to confirm the reproducibility of the observed effects. Sham stimulation was applied to rule out subjective manifestations. A new stimulation trial was initiated only after returning to the baseline pattern of the SEEG trace or after the clinical signs and symptoms had ceased. Clinical effects that patients recognized to be part of its typical ictal semiology were excluded from the analysis.

Clinical effects were then classified into six large categories according to the common descriptors of ictal behavior presented in the most recent ILAE classification (Fisher et al., 2017) to date, adapted to effects evoked by intracranial stimulation. For instance, we have not considered emotion‐related effects, as none of our patients reported them. Each large category was then subdivided into multiple levels to better characterize each clinical effect (Table 2).

TABLE 2

Three‐level categorization of the clinical effects evoked by electrical stimulation in our patients, derived from Fisher et al. (2017)

Level 1	Level 2	Level 3	Level 1	Level 2	Level 3
Cognitive	Learning	Dyscalculia	Motor	Automatisms	Laughing
	Language	Aphasia			Grimace
		Verbal fluency			Grasping
	Attention	Working memory			Rubbing
		Forced thinking			Vocalization
		Responsiveness impairment		Elementary	Dysarthria
	Memory	Déjà vu			Dystonia/tonic
		Jamais vu			Clonia
		Déjàrêvé			Weakness
		Recollection			Versive
		Memory encoding/retrieval	Sensory	Auditory	Hallucination
Autonomic	Cardiac	Bradycardia			Illusion
		Tachycardia		Gustatory	Hallucination
	Respiratory	Hyperventilation			Illusion
		Hypoventilation		Olfactory	Hallucination
	Gastrointestinal	Nausea			Illusion
		Vomiting		Somatosensory	Hallucination
	Vasomotor	Pallor			Illusion
		Flushing		Vestibular	Hallucination
		Piloerection			Illusion
		Hot–cold sensation		Visual	Hallucination
					Illusion
				Pain	Hallucination
					Illusion

For example, in the sensory domain (Level 1) a great number of clinical effects were elicited while stimulating the visual network (Level 2). Most common of these responses were simple visual hallucinations, Level 3 (e.g., patients report colored or back/white spots in the contralateral visual hemifield) and visual illusions (sensation that the image “trembles,” the image “moves”). Additionally, when stimulating the frontal lobe, we elicited eye deviation contralateral to the site of stimulation or the inability to move the eyes freely during stimulation. This type of clinical responses was classified as motor, elementary, or versive (Table 2), together with head and trunk deviation.

Similarly, in the cognitive domain (Level 1) we included attention impairment (Level 2) when patients performed worst during a digit span test (working memory deficit—Level 3), if they reported forced thinking Level 3 (an intrusive thought that appear in their mind and interfere with their ability to correctly execute a simple task) or if they show slowing in counting or reading and report they cannot concentrate or focus on a specific task without clear impairment in language or calculation ability (Table 2).

Since the first 10 patients were stimulated using a pulse frequency of 50 Hz and the remaining ones using 43.2 Hz (for minimizing the interference with power lines frequency), in order to establish the equivalence in terms of the threshold for evoking clinical effect between stimulations performed at slightly different frequencies, we have performed both 50 and 43.2 Hz stimulations in a subset of six patients. The difference in thresholds was tested against the null hypothesis using a nonparametric one‐sample Wilcoxon signed rank test.

Resolving responses from stimulation artifact

AP of biphasic pulses in a train was used to recover the EEG responses during stimulation, a method successfully used for resolve eCAP during stimulation of cochlear nerve, but not applied to date for intracranial stimulation in drug‐resistant epilepsy (Hughes et al., 2017). The basis of disambiguating the electrographic responses from the stimulation artifact is the fact that polarity of the tissue response is the same, regardless of the polarity of the stimulation pulses, whereas the artifactual components follow the AP of the stimulation pulses. The artifactual components are in a linear relationship with the stimulus, whereas the cortical responses are in a nonlinear relationship with it. While previous studies implement the artifact reduction by performing time‐domain processing of the recorded signal (e.g., Caldwell et al., 2020; Trebaul et al., 2016), we did perform a frequency‐domain analysis which allows a combined, single‐step, filtering of the response simultaneously with artifact cancelation, with a very intuitive representation in time‐frequency domain. In support of this approach in the frequency domain, we have derived the analytical equations that demonstrate that the AP stimulus does not have any components overlapping with the fundamental stimulation frequency, nor with its harmonics. Also, we have performed a numerical analysis of a simulated AP train, the responses evoked by it and a comparative study of the time‐domain versus frequency‐domain approaches.

The AP stimulus can be considered as train of biphasic square pulses () that is modulated with a square waveform having half the frequency of the regular pulse frequency (), having an amplitude ±1 (Figure 2a). A single period of the rectangular biphasic pulse, centered in the origin can be written as: where is the pulse duration (per phase) and is the period of the rectangular biphasic waveform. The way it is defined, the function is an odd function.

A single period of the modulating rectangular function having half the frequency of the stimulation waveform can be written as:

It can be noted that is an even function.

The waveforms' angular frequencies are and

The biphasic and modulating waveforms can be considered an infinite repetition of the one‐period waveform: Furthermore, the Fourier transform of a periodic function can be written as (Badrieh, 2018): where are the Fourier coefficients:

For our biphasic pulse: All coefficients are imaginary, which is a consequence of the fact that the function is odd. Their absolute value is shown in Figure 2b (black lines).

Similarly, the Fourier transform of the modulating waveform is: The Fourier coefficients are: Since , , we obtain: where coefficients are null for even values of , when both sine terms are zero, as shown in Figure 2b (gray lines).

The Fourier transform of the modulated waveform can be derived using the convolution theorem:

But by using the sifting theorem of the Dirac function (Marks II, 2009), the previous equation can be written as: Since is null for even values of , the double sum describing the frequency spectrum of the stimulation waveform does not contain any terms that correspond to an integer multiple of at which the pulses are applied. By taking , Equation (13) can be written as: that corresponds to the Fourier transform of a periodic function having the angular frequency and the Fourier coefficients: which are null for even values of , since all are null for even indexes. Therefore, the spectrum of the AP waveform will only contain non‐null terms at , , , and so on. The Fourier spectrum of the AP signal is shown in Figure 2c.

By contrast, evoked responses that are normally nonlinear and have the same polarity for each stimulation pulse, will have nonoverlapping components at integer multiples of .

These theoretical considerations are confirmed by a numerical analysis of simulated signals. The results of a numerical analysis using Fast Fourier Transform (FFT) of a discrete‐time AP waveform (Figure 3a) compared to a rectified instance of it, that emulates the constant polarity of the responses, regardless of the stimulation pulse polarity, are shown in Figure 3b. In this figure, we see that the spectral lines of the rectified signal/response (green) are distinct from the spectral lines of the AP stimulation artifact (red). This behavior is supported by the analytical equations we have derived above. It can be also more intuitively explained by the modulation theorem in Fourier analysis (Papoulis, 1962). If we take into account only the fundamental frequency of the square modulating function, the modulation theorem can be written as: where is the Fourier transform of the constant polarity stimulus function

FIGURE 3

Numerical analysis of simulated signals and a comparison between time‐domain (TDA) and frequency‐domain analysis (FDA). (a) alternating polarity biphasic stimulation train; (b) calculated frequency spectrum of the stimulation train (red) and of a rectified derivative of it (green); (c) illustration of how the modulation theorem applied to first harmonic of signals explains the splitting of the spectral lines of the artifactual components, opening up a window for detecting the underlying tissue response; (d) simulated signals obtained by differentiating Gaussians having an amplitude of 1/10 of the stimulus artifact, in the range , ] and latencies in the range , ], where , ms; the shaded area indicates the blanking interval containing the residual signal artifact excluded from the TDA analysis; (e) time‐frequency decomposition of a 5 s stimulation train based on the dark‐green waveform in (d) plus the artifact; (f) time‐frequency decomposition of a train where only the artifact is present; (g) comparison of the TDA and FDA performance, for various values of a gaussian having a latency of 12.7 ms; (h) comparison of the TDA and FDA performance for various latencies of a Gaussian having a ms; (i) scatterplot of FDA versus TDA responses for the points in G and H

The effect of such modulation of a train having the angular frequency with a waveform having a frequency is the appearance of side frequencies at and (Figure 3d), and, most notably, the disappearance of the original frequency (Figure 3c, red lines). As the neural tissue response is nonlinear, that is, has the same polarity on every pulse with a repetition rate of , the spectral content of the physiological responses will have a line at the fundamental frequency (Figure 3c, green line), plus higher‐order harmonics. This is very similar to the responses to pattern reversal stimuli in SSVEP studies (Norcia et al., 2015). Given the fact that the stimulation and modulating waveforms are square, higher‐order harmonics will be present in the spectrum, but a similar argumentation can be applied to all harmonics of the signals. Our analytical equations and numerical analysis show that they will never overlap with the response's spectrum (Figure 3b). In summary, modulating the stimulation pulse polarity splits the spectral lines of the artifactual components, opening up a spectral window exposing the underlying tissue response.

To quantify the magnitude of the responses to HFS, we have therefore implemented a bandpass comb filter that extracts from the signal the frequencies corresponding to the stimulation frequency and up to fifth order harmonics, using an ideal filter having infinite side‐lobe attenuation, as implemented in Matlab (Natick, MA). The pass band of the filter is 4 Hz at each frequency. The RMS of the filtered signal was considered to be encoding the response magnitude.

We have validated this approach by comparing the results of a standard time‐domain analysis (TDA) and our frequency‐domain analysis (FDA) on simulated signals. A response was simulated by differentiating a Gaussian with specified and latency with respect to the stimulation pulse (Figure 3d). For time‐domain analysis, the RMS of the signal obtained by averaging responses to positive and negative polarity pulses is calculated. The residual stimulus artifact was excluded from TDA analysis by blanking the signal for a duration of 5 ms covering the stimulation (Alvarez et al., 2007; Hashimoto, Elder, & Vitek, 2002; Weiss, Flesher, Franklin, Collinger, & Gaunt, 2019), duration justified by the typical appearance of the stimulus artifacts in our recordings, shown later in Figure 4d. A total of 25 different signal sets were obtained by taking 5 different values in the [2.3, 11.6] ms interval and 5 different latencies in the [6.9, 18.5] ms interval. The time‐frequency map of a 5 s simulated train containing both artifact and response in a ratio of 10 (for peak amplitudes) is shown in Figure 3e, while a map of the artifact alone is shown in Figure 3f. The Figure 3e,f show spectral lines at frequencies in agreement with our analytical equations and the numerical analysis in Figure 3b. The comparison between performance of FDA versus TDA as a function of Gaussian width for a latency of 12.7 ms is presented in Figure 3g, whereas the results of a similar analysis for different latencies of a signal having ms is shown in Figure 3h. A scatterplot of FDA versus TDA for all 25 sets is shown in Figure 3i. Figures 3g–i show that FDA constantly provides a higher, constant, value of the response since it takes into account the whole signal, without removing any part of it, whereas TDA provides lower values as it involves blanking of the signal during the stimulation and shortly afterward, that may contain part of the response waveform.

FIGURE 4

(a) Recorded EEG traces (n = 112) using a bipolar montage when stimulating the pair X'11‐X'12 located in the left Rolandic cortex of patient 12, with the baseline (green) and cortical responses to direct electrical stimulation (red) obtained by comb‐filtering the raw signals (gray), as described in 2.6; the filtered signals are represented with a magnification factor of 10 compared to raw signal; (b) magnitude of the responses on the 112 recorded channels; (c) 3D view of all SEEG electrode contacts locations (12 electrodes, 168 contacts) in patient 12, with the 128 recording contacts highlighted in green; (d) time‐domain analysis of waveforms by averaging responses to opposite pulse polarities (red, blue), illustrating incomplete artifact cancelation (green); (e) time‐amplitude map of the interpulse responses, illustrating fast onset and sustained evoked responses; (f) time‐frequency map of the signal recorded on pair M'05–M'06, illustrating spectral properties in agreement with simulation in Figure 3b; (g) same as (f), but for X'01‐X'02 pair where only artifactual components are visible; (h) same as (f), but for N′05–N′06 pair, where a 50 Hz line noise component is present, distinct from the 43.2 Hz stimulation frequency

Effective connectivity associated with clinical effects

As metric for the effective connectivity we have used the in‐band responses ( and harmonics), time‐locked with AP electrical stimulation, as described in Section 2.6. The responses were analyzed over an interval extending 2 s starting 0.25 s after the stimulation onset and were baseline‐subtracted based on a matching 2 s interval ending 0.25 s before stimulation onset. To quantify the robustness of the responses, we have further subdivided the response and baseline intervals in a sequence of 0.2 s subintervals (Figure 4a) with 50% overlap. The statistical significance of the responses was assessed using a Wilcoxon rank‐sum test comparing the set of response values on the 19 subintervals over the 2 s stimulation and baseline epochs. Responses were considered significant when the Z‐score of the differences between the two sets of subintervals exceeded 3 and the Wilcoxon rank‐sum p <.05. At the patient level, connectivity matrices including responses to all stimulations evoking a particular effect were calculated, and based on them, circular diagrams and 3D connectomes of the upper quartile (Q3) of the significant responses were represented. The group level connectivity analysis, instead of matrices and circular diagrams, makes use of the more suggestive projections of the average connectivity/responses to inflated brain surfaces, that provides visual cues about spatial locations of the more connected/activated brain areas.

In the context of our study, where we would like to capture the spatial extent of the outbound effective connectivity during HFS, among many available network metrics (Rubinov & Sporns, 2010), we have chosen to calculate the outdegree for each stimulation, as the number of significantly activated sites (Z >3, p <.05). To perform a group‐level analysis, we are required to compare networks having a variable number of nodes, represented by the number of valid channels recorded in each patient (having 64 or 128 contacts implanted), after excluding the epileptogenic and artifacted channels, which are actually variable for each stimulation. We have therefore introduced a normalized outdegree (NOD), represented by the standard outdegree, normalized by the number of valid channels, and averaged NOD over all stimulations evoking a particular effect.

Intersubject averaging, parcellation, and cortical activation maps

We have followed the methodology we previously described in Oane et al. (2020). Briefly, at the patient level, cortical surface reconstructions were performed using FreeSurfer software (Fischl, 2012) and individual contacts (or pairs) were projected to the nearest vertex of the pial surface. To perform intersubject averaging, pial surface points projected on a spherical surface model are pooled across patients and projected back to the FreeSurfer's fsaverage template (Fischl et al., 1999). Finally, a fine‐granularity multi‐modal HCP‐MMP parcellation (Glasser et al., 2016) projected onto fsaverage (available at https://figshare.com/articles/HCP‐MMP1_0_projected_on_fsaverage/3498446/2) was used to label each contact location. In addition, labels corresponding to the contacts located in two important subcortical structures that were sampled in our patients, amygdala and hippocampus were manually added and handled through a different workflow not involving their surface‐projected coordinates. For visualization purposes, generic elliptical parcels were defined on the FreeSurfer medial wall.

Current thresholds maps for evoking a clinical effect were obtained by averaging stimulation thresholds of sites having the same label and representing the values on the inflated fsaverage template. Similarly, activation maps were obtained by averaging the responses at recording sites located in each parcel, for all stimulations evoking a particular effect. As inclusion criteria for creating per‐symptom activation maps we have used a minimum number of five stimulations in two different patients.

RESULTS

In the 26 patients included in our study (17 male, 9 female), we have implanted a total of 337 electrodes (13.0 ± 3.0, mean ± SD) having 4,351 contacts (167.3 ± 42.1). The location of the electrodes, AC‐PC aligned, is shown in Figure S1.

Single stimulation example

An example of the responses to a single HFS stimulation of a site (X'11–X'12) located in the left Rolandic cortex of patient 12, evoking elementary motor effects (tonic contraction of the right hand), is presented in Figure 4. Strongest connections are local, with the electrodes M' and N' covering supplementary motor area (SMA), premotor cortex (PMC), middle and posterior cingulate cortex (MCC/PCC), primary motor and sensory cortex (R and S), and superior parietal lobule (SPL), as it can be seen in Figure 4b,c. The canonical time‐domain analysis of the signal on M'05–M'06 pair, shown in Figure 4d, performed by averaging responses to positive and negative polarity pulses show a response having same polarity and a peak at about 7 ms after the start of the stimulation pulse. The time‐amplitude map of the interpulse responses aligned with the onset of the stimulation train on the horizontal axis and each pulse on the vertical axis show that the onset of the responses is fast and the evoked potentials are steady, as it can be seen in Figure 4e. The time‐frequency plots in Figure 4f–h support the analytical equations we have derived and the numerical analysis we have performed, described in detail in Section 2.6. The time‐frequency decomposition of the signal in Figure 4d illustrates a situation where both lines corresponding to the artifactual components and the response are clearly visible and they follow the modeled distribution, including the relative amplitude of the higher‐order harmonics. The spectral lines corresponding to the artifact are at , , , and so on (Figure 4g, signal with mostly stimulus artifact), while the spectral lines corresponding to the response are located at , , , and so on (Figure 4f,h). Little out‐of‐band activation is noticed in Figure 4f–h, which might be a fact that we have used the lowest current that evokes the clinical effect ( mA). Figure 4H showing a signal affected by line noise, supports our choice for a stimulation frequency that is different from the main line noise source at 50 Hz. In addition, the entrainment of the activity shown in complex space and quantified by the phase locking value (Amengual, Vernet, Adam, & Valero‐Cabré, 2017; Lachaux, Rodriguez, Martinerie, & Varela, 1999), as illustrated in Figure S2, further supports the analysis of the steady‐state evoked potentials we have performed, through filtering the signal in a narrow frequency band around the stimulation frequency and its harmonics. The instantaneous phase vectors in Figure S2c demonstrate a quick and steady locking of the responses to the stimulation, compared to the baseline (Figure S2b).

Single patient connectivity example

The connectivity analysis for 13 stimulations evoking somatosensory effects in a single patient (P12) is presented in Figure 5. Stimulations in the K', M', P', R', S', V′, X' electrodes (3D view of the electrodes shown in Figure 4c) located in primary sensory and motor cortex, superior and inferior parietal lobule, rolandic and parietal operculum, premotor cortex, precuneus and posterior cingulate cortex elicited significant activations on 102 contact pairs located in the primary somatosensory cortex, anterior and posterior insula, parietal operculum, primary motor cortex, premotor cortex, middle and posterior cingulate cortex, inferior parietal lobule and precuneus.

FIGURE 5

Example of the connectivity analysis for 13 stimulation sites that evoked elementary somatosensory effects in patient 12. (a) Connectivity matrix. Connections not sampled are shown in dark gray; (b) circular diagram of the third quartile of the significant connections (Z‐score >3, p <.05); (c) 3D representation of the connections in (b)

Group analysis

In all patients, a total of 1840 pairs of contacts were stimulated. About one third (614 stimulations) evoked a clinical effect, while 1,226 did not. After excluding artifacted or epileptogenic contacts, on 2,613 pairs of contacts implanted in all 26 patients, 46,211 responses were recorded, of which 12,547 (21.9%) were significant (Z‐score >3, p <.05). The histogram of the number of significant responses as a function of voltage is given in Figure S3a. The mean evoked voltage was 7.6 μV (significant responses only). The amplitude of the responses decreases following a power‐like law with Euclidean distance from stimulation site (Figure S3b). The difference between current thresholds for evoking clinical effects of the 43 and 50 Hz stimulations performed at the same sites is shown in Figure S4. There was no significant difference between the two thresholds, Wilcoxon signed rank test p >.55. Additionally, we have validated the FDA approach versus TDA in all 16 patients where 43.2 Hz stimulations were performed, results being shown in Figure S5. The FDA and TDA values are highly correlated (Spearman r = .99, p <.0001). In agreement with the results of the computational analysis shown in Figure 3g,h, the FDA responses are systematically higher than the TDA responses (slope of Deming linear regression 1.16), as FDA considers the signal in its entirety, whereas TDA blanks the signal for 5 ms around the stimulation artifact.

To capture the overall spatial characteristics of the responses, in Figure 6 we have plotted the histograms of the total number of responses and significant responses as a function of the Euclidean distance from the stimulated site for all stimulations evoking effects (a,b) and per‐effect (c,d). The distribution of the number of significant activations had a mean distance of 30.4 mm, which is considerably lower than 42.0 mm for all responses (Figure 6a). Based on these results, we have defined for the subsequent analysis short range connections as having a distance from stimulation site d <30 mm and long‐range connections as having d ≥30 mm. The percentage of significant responses, independent of the spatial distribution of contacts, follows a power‐like law (Figure 6b). When analyzing per‐symptom response distributions, we see that the motor effects (L1 classification) have a shorter activation range (mean 28.8 mm), which is statistically different from the value for the other shown classes of effects (autonomic p <.05, others p <.001). Also, based on L2 classification, elementary motor (28.7 mm) has the lowest mean distances, statistically different from others (p <.05).

FIGURE 6

(a) Histogram of the number of responses recorded on all implanted contacts, and of the number of significant responses (Z‐score >3, p <.05) as a function of the Euclidean distance from the stimulated site; (b) histogram of the ratio between significant nonsignificant responses shown in a; (c,d) histogram of the significant responses as a function of the distance from the stimulated site for L1 (c) and L2 (d) effect classification. Only top 5 effects in terms of number of evoked symptoms have been represented

Stimulated sites (n = 1840) were covering most of the cortical surface (Figure 7a). Sites that elicited a clinical effect (n = 614) are primarily located in the motor and sensory (somatosensory, visual), insular and peri‐sylvian cortical areas, as it can be seen in Figure 7a. Locations that did not elicit a clinical effect up to the maximum allowable intensity (3 mA) are included in associative cortices like prefrontal cortex, parietal and temporal neocortex. Mean intensity from stimulations that elicited a clinical effect was 1.03 mA. Among the most activated regions are the motor and prefrontal cortex (Figure 7b). Areas that exhibit short‐range activations are the ones most responsive to stimulation: sensorimotor cortex, insular and opercular cortex, supplementary and premotor cortex, mesial parietal cortex, mesial temporal structures and primary visual areas (Figure 7c). On the other hand, areas associated to distant high activations are the prefrontal cortex, the hippocampus and several subregions of the temporo‐occipital cortex (Figure 7d).

FIGURE 7

Current threshold (a) for evoking clinical effects and response maps (b–d) for all 614 stimulations in 26 patients. Sites evoking clinical effects were shown using green circles (a–d), while the others were shown using empty circles (a); (b) response maps for all stimulations and all distances between stimulation site and recorded site; (c) response maps for short‐range connections (Euclidean distance from stimulated site <30 mm); and (d) for long‐range connections (≥30 mm)

Per‐symptom connectivity

By selecting the subset of stimulations that evoked specific effects, and associated responses, we were able to create per‐effect threshold and activation maps. As an example, we have shown in Figure 8 the maps for auditory effects (L2 classification, according to Table 2), evoked by 42 stimulation sites in 12 patients. These were elementary hallucinations or illusions when patients report hearing high pitched sounds, their voice in echo or that the environmental sounds are less intense. These auditory effects are most of the time bilateral, but in a few cases patients describe them to be lateralized to the ear contralateral to the side of stimulation. Mean stimulation intensity is 0.96 mA with lowest threshold in the posterior insular cortex, long gyri (Figure 8a). Maps show mostly a local activation in the posterior insula and temporal operculum with long distance activations in the frontal operculum, anterior insula and temporal‐posterior neocortex (Figure 8b–d). Low amplitude responses are also elicited in the temporo‐basal regions, anterior temporal neocortex, inferior parietal lobule, frontal and parietal operculum, mesial prefrontal and parietal cortex.

FIGURE 8

The thresholds for the 42 stimulated sites evoking auditory effects (a); All evoked activations are shown in b, whereas local (Euclidean distance from stimulated site d <30 mm) and distant (d ≥30 mm) ones are shown in c and d, respectively

The maps for all effects we have observed in our 26 patients, grouped by classification level, are presented in Figures S6–S28.

Our data does not capture the connectivity between pairs of recorded contacts, since we are only calculating the outbound effective connectivity between stimulated contacts and recorded pairs. Among the large variety of network measures (Rubinov & Sporns, 2010), the only measure that we consider to be relevant in this context is the outdegree. Per‐effect (L1 classification) mean NOD is shown in Figure 9, with motor exhibiting the largest number of outbound connections.

FIGURE 9

The mean normalized outdegree for stimulations evoking effects according to the L1 classification. The number of stimulations evoking each class of effects is shown on the bars

DISCUSSION

Functional and effective connectivity estimates using electrophysiological signals can be based on complex‐space, amplitude, or phase methods, linear, nonlinear or high‐complexity (mutual information, directed transfer function, Granger causality, dynamic causal modeling and other) connectivity metrics (see Rossini et al., 2019 for a review). Most of the studies cover resting‐state connectivity, given the difficulty in capturing and interpreting the dynamic processes taking place during the execution of tasks.

To date, task‐related functional connectivity has been studied primarily using fMRI (Glover, 2011). However, there is no experimental paradigm that can guarantee the focal activation of a cortical area, like intracranial stimulation does. Also, functional MRI connectivity measures are based on indirect evidence of the cortical activations, using BOLD signals that have a low temporal resolution. Our study uses for both stimulation and quantification of response the native electrophysiological signals the brain. Task‐based fMRI can only test the cortical activations for a finite set of experimental paradigms that the researchers can think of or that can be implemented. The inferred connectivity is often dependent on the choice of the seeds/ROIs, which is not always guided by objective or independent criteria. Our approach is based on an unbiased exploration of the cortex's function by applying standardized electrical stimulation and monitoring quasi‐instantaneous EEG responses to HFS stimulation.

While other studies present detailed maps of the cortical locations evoking effects, the clinical value of the stimulation current thresholds is not discussed. We highlight that cortical regions that elicited a clinical effect at lowest threshold intensity overlap with those area considered to be eloquent cortex: e.g., primary somatosensory area, visual cortex or language areas (Luders, Schuele, & McIntyre, 2008). The eloquent cortex overlap with the areas that evoked clinical effects at low stimulation thresholds (<1 mA), suggesting that there may be a relation between increased cortical excitability and brain function, and stressing the need of careful consideration when including these areas in the resection volume, to avoid potential deficits.

While a recent study by Perrone‐Bertolotti et al. (2020) analyzes the language networks by looking at out‐of‐band activations, specifically in the high‐frequency 70–150 Hz range (without including any harmonics of the stimulation frequency), we specifically analyze the in‐band responses at the stimulation harmonics and up to fifth harmonic, where most of the spectral power of the signal is located. This does not capture the secondary activations, but the extent of the primary brain networks effectively entrained by the stimulation, time‐locked to it, in the initial 2 s from the stimulation onset.

Key to this study is the implementation of methodology that completely disambiguates responses from the artifactual components by combining a modification of the stimulation protocol consisting in alternating the polarity of the pulses in a train with a response analysis in the frequency domain. The classical time‐domain analysis of the responses to AP HFS through response averaging (Hughes et al., 2017) results in an incomplete cancelation of the stimulation artifact, as visible in our data in Figure 3g. The frequency‐ domain approach we introduce addresses several issues of the time‐domain analysis, that includes pulse‐to‐pulse variability of the stimulation artifact due to misalignment of the stimulator's and EEG recording system's clocks (Caldwell et al., 2020; Hashimoto et al., 2002), unequal excitation between polarities and unequal latency between polarities (Hughes et al., 2017). All these factors lead to a variation of the phase of the recorded signal, a broadening of the spectral lines and a modified ratio between various harmonics, but do not alter the basic frequency content of the signal. Since, in our approach, the response magnitude is calculated based solely on the absolute value of the complex values resulting from the frequency filtering, all the effects of the phase variations are nullified.

This methodology allows highlighting clinical effects' networks, which goes beyond the traditional view according to which “brain processing was mainly thought in a localizationist framework, in which one given function was underpinned by a discrete, isolated cortical area”—Herbet and Duffau (2020). In that respect, our study bridges the concept of functional segregation, referring to the localization of brain function, with the concept of functional integration, referring to the involvement of several specialized areas to accomplish a particular function (Friston, 2011), as a necessary step in coming up with novel, highly complex, meta‐networking models of the brain (Herbet & Duffau, 2020).

Small differences between the mean spatial extent of the networks underlying different symptom classes have been observed (29.0–31.9 mm), though some of the differences were found to be statistically significant. This shows that the spatial characteristics of the activated networks are quite similar, regardless of the evoked symptoms. Notably, there is a statistically significant difference between motor and cognitive/multimodal effects (p <.005), the first ones being associated with the shortest distance connections (Figure 6c,d). On the other hand, cognitive effects are associated with the smallest percentage of NOD connections (Figure 9). We highlight here that the brain's effective connectome is organized in stronger short‐range connections, most of them related to elementary functions (e.g., motor, sensory) and an apparently weaker long‐range connections associated with complex functions (e.g., cognitive Figure S7, working memory Figure S28). This view is supported by previous studies, which point out that long‐distance connections facilitate rapid and efficient interareal communication (Betzel & Bassett, 2018). The authors also suggest that distant connections exhibit redundancies that could be the reason why associative cortices (characterized by long range connections Figure 7d) do not display clinical effects during focal electrical stimulation, and do not generate a deficit after a cortical resection. However, a case are those cortical regions where autonomic responses are elicited, that have the most widespread (Figure 9) and long‐range connections (Figure 6c). Despite the high metabolic cost, it seems that long‐range connections play critical functional roles (Betzel & Bassett, 2018) and the large number of outgoing connections could be explained by the fact that most of the human behavior is accompanied by autonomic responses; increase in heart rate while performing cognitive tasks or motor exercises (Critchley, Corfield, Chandler, Mathias, & Dolan, 2000). One has to keep in mind that the spatial distribution of the responses in Figure 6 is not uniquely determined by the brain networks' architecture but is shaped by the spatial distribution of the contacts, based solely on the surgical plans of each patient. With that respect, the ratio between significant and all recorded responses, shown in Figure 6b, is a measure that is independent of the contacts' spatial distribution, and, as expected, follows a power‐like law (Trebaul et al., 2018).

The per‐effect activation maps as presented in Figure 8 and in Figures S6–S28 constitute a prototype of an atlas of networks associated with specific brain functions. Their clinical value for the pathological brain refers to the ability to create network models of functional cortical organization which could help us understand how brain damage lead to clinical expression (Hallett et al., 2020). In this way, they can also contribute to understanding the patterns of spatial propagation of ictal activity, based on the timeline of ictal semiology. Epilepsy is increasingly considered a network disorder (Kramer & Cash, 2012), therefore understanding pathological and physiological networks, resting state or associated with brain function, is equally important in understanding this disease.

The limitations of this study refer primarily to the sparse and nonuniform sampling of the brain with intracranial electrodes. There may be areas that are activated by electrical stimulation, but we fail to detect their activation since they are not covered with electrodes. With that respect, the results of the study should be cautiously regarded as a sampling of the brain networks underlying clinical effects, not necessarily as a method proving what these networks are. The number of stimulations evoking each class of effects is highly variable in our subjects, as indicated in Figure 9, which in combination with patient's sampling bias, might possibly be standing behind some unexpected results. A very important fact to take into consideration is that for depth electrodes like the ones we used, the activation of the white matter tracts passing by the electrodes, part of different networks, cannot be ruled out (Duffau, 2015). Electrical stimulation propagates through the underlying anatomical connections of the stimulation and response sites (Conner et al., 2011; Donos et al., 2016; Silverstein et al., 2020) therefore it may be possible that the stimulation‐evoked responses extend beyond the effect network. In other words, the effect network is included in the activation map, but the activation map might include additional brain structures, which are not necessarily related to the clinical effect. While it might have been of interest to study the relationship between activations and anatomical connections between stimulated and response areas (Conner et al., 2011; Donos, Mălîia, et al., 2016; Silverstein et al., 2020), since high‐quality diffusion imaging was not available in all patients, we believe that an atlas‐based approach would not have provided robust results, particularly for our patient sample size and fine‐granularity of the HCP‐MMP parcellation we have used. Our study aims to primarily introduce the methodological approach for highlighting the networks activated by DES, and introduce the potential clinical value of the activation maps, and not to elucidate all neurophysiological mechanisms underlying the activations.

One has to note that clinical effects, particularly nonmotor ones, can be complex and represent a subjective phenomenon experienced only by the patient. In these cases, no objective manifestation is present, therefore collection and classification of these results are based on patients' ability to explain what they feel, which is an inherently subjective process.

Perhaps one potential weakness of the study is the fact it based on a single‐trial paradigm and that the statistical test for assessing the significance of the responses might be considered weak, being based on a set of 19 responses on 0.2 s sub‐intervals for each condition (response/baseline). This is dictated by the maximal duration of the stimulation of 5 s and the fact that some stimulations are terminated prematurely when ample symptoms are evoked, limiting the practical duration of the analysis interval to 2 s. One other limitation refers to the fact that although the frequency‐domain analysis of the responses to AP stimulation pulses is quite robust in theory, there may be some practical aspects of the signal collection that may result in a broadening of the spectral lines, such that the tails of the spectral peaks of the artifactual components may extend into the response frequencies of interest. Such effects can be the result of different factors, like variation of the phase or shape of the artifactual components due to the mismatch and drift in the clock of the stimulator and of the EEG recording system (Caldwell et al., 2020; Hashimoto et al., 2002), or nonlinearities in the ADC conversion, particularly when signal level is approaching maximal range of the input. The high sampling rate (4,096 Hz) and the quality of the recording systems used in our study is supposed to minimize these issues. These spurious contributions to the response, though small in amplitude, may be resulting in false positive detections in the terms of response significance, as they are sustained over the 19 sub‐intervals covering the response intervals. As concerns the method for visualizing the thresholds and responses over the inflated brain, one must keep in mind that the process of cortical surface reconstruction, mapping the contacts to the surface and group analysis may be affected by errors (Fischl, 2012), resulting in an approximate localization and consequently parcellation assignment of the stimulated or recorded contacts over the brain template.

CONCLUSIONS

Key to enabling visualization of the networks associated with clinical effect is a methodology that by modulating the polarity of the pulses in a stimulation train opens up a spectral window allowing the analysis of the responses during intracranial high‐frequency stimulation. This allowed us to characterize the organization and spatial extent (local or distant) of the clinical effects' networks. In addition, our study provides whole‐brain stimulation threshold maps for evoking clinical effects, highlighting the resilience of the associative cortex to stimulation and responsiveness of the granular cortex.

CONFLICT OF INTEREST

Andrei Barborica is also Vice‐President and Chief Technological Officer of FHC Inc., the manufacturer of the electrical stimulator and stereotactic fixture used in this study. The other authors have nothing to disclose in relation to this work.

AUTHOR CONTRIBUTIONS

Andrei Barborica: Conceptualization, methodology, validation, software, formal analysis, writing – original draft, visualization, supervision, funding acquisition; Irina Oane: Conceptualization, validation, investigation, writing – original draft; Cristian Donos: Writing‐review and editing; Andrei Daneasa: Investigation; Felicia Mihai: Data curation, formal analysis; Constantin Pistol: Data curation; Aurelia Dabu: Investigation; Adina Roceanu: Resources; Ioana Mindruta: Conceptualization, validation, investigation, resources, writing review and editing, supervision, funding acquisition.

ETHICS STATEMENT

The study has been performed under Emergency University Hospital Bucharest ethical committee approval 32483/27.06.2018 and Bucharest University ethical committee approval CEC 45/11.06.2020.

Patient consent

All patients signed a written informed consent, in accordance with the Declaration of Helsinki, for the recordings, stimulations and data sharing procedures.

Appendix S1: Supporting Information

Click here for additional data file.