Simultaneous human intracerebral stimulation and HD-EEG, ground-truth for source localization methods

Precisely localizing the sources of brain activity as recorded by EEG is a fundamental procedure and a major challenge for both research and clinical practice. Even though many methods and algorithms have been proposed, their relative advantages and limitations are still not well established. Moreover, these methods involve tuning multiple parameters, for which no principled way of selection exists yet. These uncertainties are emphasized due to the lack of ground-truth for their validation and testing. Here we present the Localize-MI dataset, which constitutes the first open dataset that comprises EEG recorded electrical activity originating from precisely known locations inside the brain of living humans. High-density EEG was recorded as single-pulse biphasic currents were delivered at intensities ranging from 0.1 to 5 mA through stereotactically implanted electrodes in diverse brain regions during pre-surgical evaluation of patients with drug-resistant epilepsy. The uses of this dataset range from the estimation of in vivo tissue conductivity to the development, validation and testing of forward and inverse solution methods.

2 Scientific Data | (2020) 7:127 | https://doi.org/10.1038/s41597-020-0467-x www.nature.com/scientificdata www.nature.com/scientificdata/ robustness 6,7 . That is, simulated electrical activity is placed inside a realistic volume-conductor model and projected onto the scalp surface in order to be used as input data for source localization algorithms, which are then tested on their ability to reconstruct the origins of these signals. Another common methodology is to try localizing functional activity whose origins are inferred from other imaging modalities 8 (i.e. fMRI during somatosensory stimulation). However, simulations lack realism and cross-modal functional mapping lacks spatial precision and can introduce relative biases in spatial arrangement due to the different nature of the signals.
A fundamental element to fill this gap could be offered by stereo-electroencephalography (sEEG), obtained from drug-resistant epileptic patients using stereotactically implanted electrodes. Once surgically implanted, patients are monitored continuously for several days to have one or more seizures recorded. During this time, sessions of intracortical stimulation are performed in order to induce habitual seizures and to provide a map of the physiological functions of the implanted sites [9][10][11][12][13][14] . This procedure implies that a brief current pulse is injected between two adjacent leads, producing an electrical artifact whose localization can be accurately determined. When combined with simultaneous scalp EEG, this procedure is capable of generating real data of scalp recorded electrical signals originating from precisely known locations inside the human brain, and thus represents an ideal benchmarking scenario for validating and comparing both forward and inverse solution methods.
In line with this, the aim of this paper is to provide a consistent dataset of high-density scalp EEG recordings performed during the stimulation of intracortical leads. It contains the anonymized MRIs necessary to build forward models, the surfaces and forward models created using the subjects' original MRIs, the spatial and anatomical information of the stimulated sites, and EEG data from 256 channels with digitized positions. As a further element, stimulations were performed at different current intensities, so as to favor not only a comparative performance across different topographical regions, but also an estimation of the role that the intensity of a source activity plays in its localization accuracy. The value of this dataset is also increased by the dense sampling of the scalp, which allows spatial down-sampling procedures to test the performance of inverse solution algorithms under a montage-dependent perspective.
In order to demonstrate the validity and wide range of possible uses of this dataset, we performed six different analysis. First, we tested the performance of three widely used inverse solution methods, employing various montages and parameters' configurations, and tested the best reachable performance. Second, we examined how misselection of parameters affected localization accuracy. Third, we analyzed the spatial dispersion of the computed solutions across methods and montages. Fourth, we assesed the spatial profile of the observed localization errors. Fifth, we characterized the relationship between localization errors and depth of the stimulated sites. Finally, we evaluated how different MRI anonymization procedures influence source localization results.
To the best of our knowledge, Localize-MI would be the first dataset providing the neuroscientific and technical community with ground truth to validate the efficacy of forward and inverse solutions on EEG data, and to systematically evaluate the factors mostly contributing to the overall process accuracy.

Methods
Participants. Seven subjects (F = 4) participated in the study (X age = 35.1; sd age = 5.4). A total of 61 sessions were obtained (X sessions per subject = 8.71; sd sessions per subject = 2.65). All subjects were patients undergoing intracranial monitoring for pre-surgical evaluation of drug-resistant epilepsy ( Table 1). All of them provided their Informed Consent before participating, the study was approved by the local Ethical Committee (protocol number: 463-092018, Niguarda Hospital, Milan, Italy) and it was carried out in accordance with the Declaration of Helsinki. All subjects underwent surgical or thermocoagulation procedures with less than two years of follow-up time, therefore proper Engel surgical outcome classification scores 15 were not available. However, their corresponding values would be Ia for all of them; with the exception of sub-02, where it would be IIa, and sub-05, where the procedure was carried out with less than 2 months of follow-up time and we cannot provide an approximative score.
Electrical stimulation. Intracranial shafts were implanted using a robotic assistant (Neuromate; Renishaw Mayfield SA), with a workflow detailed elsewhere 13 . The position of the implanted electrodes was decided exclusively following clinical needs. Stimulation sites were chosen in collaboration with the epileptologist in charge of the patients. We selected contacts that were located in anatomically normal brain regions, outside the epileptogenic zone and without pathological sEEG activity. Electrical currents were delivered through platinum-iridium semiflexible multi-contact intracerebral electrodes (diameter: 0.8 mm; contact length: 2 mm, inter-contact distance: 1.5 mm; Dixi Medical, Besançon, France). Single-pulse biphasic currents lasting 0.5 ms were delivered at intensities ranging from 0.1 to 5 mA (number of sessions: 0.1 mA = 22; 0.3 mA = 17, 0.5 mA = 8; 1 mA = 9; 5 mA = 5) through pairs of adjacent contacts by a Nihon-Kohden Neurofax-100 system (Fig. 1). The stimulation frequency (i.e. number of pulses per second) was of 0.5 Hz when stimulating at 1 and 5 mA and 1 Hz otherwise (with the exception of 3 sessions at 1 mA on which the stimulation frequency was 1 Hz). A total of 60 trials were obtained from each stimulation site when stimulating at 0.1, 0.3 and 0.5 mA, and a total of 40 when stimulating at 1 and 5 mA (Fig. 2). We chose to use the stimulation artifacts instead of specific brain rhythms because of the precise spatial location of the former, given that in the case of brain rhythms, their generators might not be exactly at the location of the intracranial contact, and may therefore bias the estimation of the methods' accuracy. EEG recordings. EEG Recordings were performed at the end of the pre-surgical evaluation. The EEG cap was sterilized and, after the protective bandage was removed, the skin was disinfected and the cap placed in position. At the end of the recording session the skin was disinfected again. The whole procedure was carried out by neurosurgeons using sterile technique. EEG signals were recorded from 256 channels (Geodesic Sensor Net; HydroCel CleanLeads) sampled at 8000 Hz with an EGI NA-400 amplifier (Electrical Geodesics, Inc; Oregon, USA), using a custom-built acquisition software written in C++ and Matlab (The Mathworks Inc.), based on www.nature.com/scientificdata www.nature.com/scientificdata/ EGI's AmpServerPro SDK. All software filters were disabled during acquisition. The spatial locations of EEG electrodes and anatomical fiducials were digitized with a SofTaxicOptic system (EMS s.r.l., Bologna, Italy), coregistered with a pre-implant MRI (Achieva 1.5 T, Philips Healthcare).

Electrode localization. The location of the intracranial electrodes was assessed registering the post-implant
CT (O-arm 1000 system, Medtronic) to the pre-implant MRI by means of the FLIRT software tool 16 . The position of every single lead was assessed with respect to the MRI using Freesurfer 17 , 3D Slicer 18 and SEEG assistant 19 . When the pre-implant MRI and the EEG digitization MRI were not the same, contacts positions were transformed from the SEEG space to the EEG space using an affine transformation between MRIs calculated employing the ANTs software 20 . Normalized contacts' coordinates were estimated by performing a non-linear registration between the subject's skull stripped MRI and the skull-stripped MNI152 template 21 (ICBM 2009a Nonlinear Symmetric) using ANTs' SyN algorithm. Contact positions were plotted on a flatmap of the MNI152 template built using Pycortex 22 , by projecting each contact's coordinates to the closest vertex of the brain surface reconstruction. The accuracy of the normalization procedure was verified by visual inspection.
Data preprocessing. Raw data were imported and preprocessed in Python employing custom-built scripts and the MNE software 23,24 . Continuous data were high-pass filtered at 0.1 Hz (FIR filter; zero phase; Hamming window; automatic selection of length and bandwidth). Data from two subjects (sub-05 and sub-07) were also notch filtered at 50, 100, 150 and 200 Hz (FIR filter; zero phase; Hamming window; bandwidth = 0.1 and automatic length selection) due to considerable line noise. Bad channels were identified by visual inspection (i.e. flat channels, presence of artifacts, etc.). Next, epochs were generated from −300 ms to 50 ms with respect to the stimulation electrical artifact and baseline corrected (mean subtraction method, from −300 ms to −50 ms). The baseline period was specifically chosen to avoid any possible contamination by cortico-cortical evoked responses from previous trials, even with the fastest stimulation frequency 25 . Bad epochs were identified by visual inspection and rejected. Given that EGI's trigger channel is sampled at 1000 Hz, which introduced jitter between the onset of the trigger and the onset of the stimulation, epochs were fine-aligned by matching the peaks of the stimulation artifacts within sessions. All good epochs were saved in MNE's fif format in the interval between −250 and 10 ms and subsequently converted to BIDS format 26,27 . Source localization. The source localization procedure was carried out using the MNE software. Surface reconstructions were obtained with Freesurfer and a 3-layer Boundary Element Method (BEM) model was created with 5120 triangles and conductivities set to 0.3, 0.006 and 0.3 S/m, for the brain, skull and scalp compartments respectively. Source spaces were created with 4098 sources per hemisphere. Epochs were re-referenced to the average of all good channels and covariance was estimated with automated method selection 28 . Subsequently, epochs were averaged and cropped from −2 to 2 ms with respect to the stimulation artifact. Inverse solutions were calculated with three different methods: Minimum Norm Estimate (MNE), dynamic Statistical Parametric Maps (dSPM) and exact Low Resolution Electromagnetic Tomography (eLORETA) 5,29-31 . These methods were selected based on two criteria. The chosen methods had to be (1) widely used by researchers, in order to be representative of currently used methods, and (2) implemented in an open-source and free to access software platform, in order to favor its access and facilitate reproducibility.
Various parameter configurations were assessed. The regularization parameter was set as 1/SNR 2 with SNR set to 1, 2, 3, and 4. The depth and loose weighting parameters varied between 0.1 and 1 in 0.1 steps. Four different EEG montages were tested: all good channels, and channels corresponding to EGI's 128, 64 and 32 montages. www.nature.com/scientificdata www.nature.com/scientificdata/ When a channel selected for the subsampled montage was marked as bad, we replaced it by its closest neighbour. A total of 4800 solutions were calculated for each session.
The Euclidean distance between the coordinates of the center of the pair of stimulating contacts and the coordinates of the maximal activation in the source estimates were computed as well as the distance on each spatial axis (left-right, anterior-posterior and inferior-superior) as measures of accuracy. We then computed the best solution across all montages and parameter's configurations.
We also calculated number of sessions on which each method and montage reached the minimum distance and the proportion of solutions for each of these sessions on which they were able to reach it (i.e. the number of solutions for a session and method or montage on which it reached the minimum distance divided by the total number of solutions computed for that montage or method).
In order to further evaluate the performance of each method and montage we also computed the Spatial Dispersion metric 32,33 , which is calculated as shown in Eq. 1:  www.nature.com/scientificdata www.nature.com/scientificdata/ Finally, in order to evaluate the performance of each method as a function of the depth of the stimulated site, we performed a mixed-effects linear regression analysis with distance to the stimulation site as dependent variable, distance to the skin as predictor and subject as random factor (intercept). The mixed-effects approach was chosen due to the nested nature of the data (i.e. stimulation sites within subjects) 34 . For one method (eLORETA) the mixed-effects model resulted in a "singular fit" due to a lack of variance in the random intercept and we therefore performed a standard linear regression analysis. The distance to the skin was calculated as the minimum distance from the position of the stimulated contacts and the skin surface obtained with the watershed algorithm used for the BEM model. Marginal R 2 was used to calculate the variance explained by the models 35 , and Adjusted R 2 for the standard linear regression analysis. MRI anonymization. MRIs were anonymized employing two different tools: Pydeface (https://github.com/ poldracklab/pydeface) and Maskface 36 (Fig. 2e). In order to investigate the influence on source localization results of the geometrical distortions induced by the anonymization procedures, we recreated the forward-models with the anonymized MRIs and computed the inverse solutions of all the parameters' configurations that reached the minimum distance of each session. We then compared the distances to the stimulation sites obtained with the anonymized MRIs with the ones obtained with the original ones.

Data Records
The Localize-MI dataset is available at the Human Brain Project platform 37 (https://doi.org/10.25493/NXN2-05W) and at G-Node 38 (https://doi.org/10.12751/g-node.1cc1ae). The dataset comprises high density-EEG data from a total of 61 sessions, obtained from 7 subjects (Online-only Table 1). In addition, it includes the spatial locations of the stimulating contacts in native MRI-space, MNI152-space and Freesurfer's surface-space, and the digitized positions of the 256 scalp EEG electrodes. It also contains the surfaces used for creating the BEM models, the pial and inflated surface reconstructions created with the subjects' original MRIs, as well as the source-spaces and forward-models from them derived. Furthermore, it includes the anonymized MRIs of each subject. www.nature.com/scientificdata www.nature.com/scientificdata/ technical Validation Methods, montages and parameters. The minimum distance between the stimulation sites and the location of the maximum current values was between ~2 and ~12 mm when optimal parameters were selected (X minimum distance = 5.39 mm; sd minimum distance = 2.61, min minimum distance = 2.30, max minimum distance = 12.16; Fig. 3a). Instead, when all parameters' configurations were considered, the distance between the stimulation site and the location of the maximum current values was generally between ~2 mm and ~50 mm (Fig. 3b,e).
The mean of the proportion of solutions for each session on which each method reached the optimal solution was 0.02 for MNE on a total of 11 sessions, 0.02 for dSPM on a total of 30 sessions and 0.06 for eLORETA on a total of 32 sessions (Fig. 3c). The mean Spatial Dispersion was 38.6 for MNE, 40.7 for dSPM and 40.6 for eLO-RETA (Fig. 3d).
The mean of the proportion of solutions for each session on which each montage reached the optimal solution was of 0.06 for all channels on a total of 15 sessions, 0.05 for 128 channels on a total of 22 sessions, 0.04 for 64 channels on a total of 26 sessions, and 0.04 for 32 channels on a total of 27 sessions (Fig. 3f). The mean Spatial Dispersion was 38.6 for all channels, 38.4 for 128 channels, 40.8 for 64 channels and 43.5 for 32 channels (Fig. 3g).
The differences between the stimulation site and the location of the maximum current value of the solutions that reached the best solution for each session were approximately centered around zero and symmetrical across the three spatial axes (Fig. 3h). www.nature.com/scientificdata www.nature.com/scientificdata/ Finally, the mixed-effect linear regression analysis (Fig. 3i) showed that the performance of MNE was remarkably related to the depth of the stimulation site, with positive slope and a high coefficient of determination (β = 0.70, R 2 = 0.71). Conversely, for dSPM the relationship was negative, with a relatively low coefficient of determination (β = −0.27, R 2 = 0.27). Finally, for eLORETA the slope was positive, but the coefficient of determination was low (β = 0.16, R 2 = 0.18).

MRI anonymization.
The distance between the stimulation sites and the location of the maximum current values remained equal in a relatively large number of solutions when employing the anonymized MRIs for the calculation of the forward models, with both anonymization methods (% equal deface = 0.89; % equal maskface = 0.89). However, a number of them proved to produce different results.

Usage Notes
The Localize-MI dataset is provided in BIDS format and contains all the necessary information to allow researchers to perform their analysis on any software. However, please note that, at the time of publication of this article, the BIDS specification for Common Electrophysiological Derivatives has not been established yet and therefore the dataset structure might not be compatible out-of-the-box with all software. However, adjusting the structure for specific purposes should be straight-forward and, importantly, once the specification will be published, we will update the database in order to conform to it. Interactive scripts of usage demonstration are provided as part of the repository accompanying this article.
This dataset has multiple potential uses, for instance: estimating in-vivo tissue conductivities; evaluating the impact of different forward-models on inverse solutions; developing, validating and testing different inverse solution methods; studying interactions between forward and inverse solution methods; performing linear combinations of stimulation sessions in order to test the ability of diverse methods to retrieve the correct sources; etc.
It is worth mentioning that the artifacts generated by intracranial stimulation are non-physiological, therefore generalization of results to physiological signals should be done conscientiously. Also, in some cases, the tails of the intracranial shafts, which protruded from the scalp, precluded the contact with the skin of a number of EEG electrodes. Nevertheless, the analysis performed revealed good localization accuracy, demonstrating that this was not an issue. Another limitation corresponds to the fact that anatomical areas sampled tend to be clustered within subjects, which should be taken into consideration when performing topographical analysis. However, the Localize-MI dataset will be extended with data from new subjects in the future, which will provide a more comprehensive spatial coverage and allow more detailed spatial analyses.

Code availability
Usage demonstration scripts and the code used for the preparation, pre-processing and technical validation of the Localize-MI dataset are publicly available at https://github.com/iTCf/mikulan_et_al_2020.