# An open-source toolbox for Multi-patient Intracranial EEG Analysis (MIA)

A.-Sophie Dubarry, Catherine Liégeois-Chauvel, Agnès Trébuchon, Christian Bénar, F-Xavier Alario

► To cite this version:

A.-Sophie Dubarry, Catherine Liégeois-Chauvel, Agnès Trébuchon, Christian Bénar, F-Xavier Alario. An open-source toolbox for Multi-patient Intracranial EEG Analysis (MIA). NeuroImage, 2022, 257, pp.119251. 10.1016/j.neuroimage.2022.119251. hal-03713056

HAL Id: hal-03713056

https://hal.science/hal-03713056

Submitted on 4 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An open-source toolbox for Multi-patient Intracranial EEG Analysis (MIA)

A.-Sophie Dubarry¹,², Catherine Liégeois-Chauvel³,⁴, Agnès Trébuchon⁴,⁵, Christian Bénar⁴*, F.-Xavier Alario²,³*

1. Aix Marseille Univ, CNRS, LPL, Aix-en-Provence, France

2. Aix Marseille Univ, CNRS, LPC, Aix-en-Provence, France

3. Cortical Systems Laboratory, University of Pittsburgh, Pennsylvania, USA

4. Aix Marseille Univ, INSERM, INS, Inst Neurosci Syst, Marseille, France

5. APHM, Hôpital la Timone, Service Épileptologie et Rythmologie Cérébrale, Marseille, France

**Corresponding author:**

A.-Sophie Dubarry (anne-sophie.dubarry@univ-amu.fr)

Aix Marseille Univ, CNRS, LPL, Aix-en-Provence, France

Laboratoire Parole et Langage. 5 avenue Pasteur, 13100 Aix-en-Provence, France

**Study funding:**

This work has been funded by the European Research Council FP7/2007-2013 Grant 263575, as well as Agence Nationale de la Recherche Grant Multimodel 2010 BLAN 0309 01. This research, carried out within the Institut Convergence ILCB (ANR-16-CONV-0002), has benefited from support from the French government, managed by the French National Agency for Research (ANR) and the Excellence Initiative of Aix-Marseille University (A*MIDEX)

**Highlights:**

* Software package for the analysis and visualization of intracranial signals

* Consistent and user-friendly graphical interface

* Scriptable package for reproducible analysis
**Abstract** (183 words)

Intracranial EEG (iEEG) performed during the pre-surgical evaluation of refractory epilepsy provides a great opportunity to investigate the neurophysiology of human cognitive functions with exceptional spatial and temporal precisions. A difficulty of the iEEG approach for cognitive neuroscience, however, is the potential variability across patients in the anatomical location of implantations and in the functional responses therein recorded. In this context, we designed, implemented, and tested a user-friendly and efficient open-source toolbox for Multi-Patient Intracranial data Analysis (MIA), which can be used as standalone program or as a Brainstorm plugin. MIA helps analyzing event related iEEG signals while following good scientific practice recommendations, such as building reproducible analysis pipelines and applying robust statistics. The signals can be analyzed in the temporal and time-frequency domains, and the similarity of time courses across patients or contacts can be assessed within anatomical regions. MIA allows visualizing all these results in a variety of formats at every step of the analysis. Here, we present the toolbox architecture and illustrate the different steps and features of the analysis pipeline using a group dataset collected during a language task.

**Keywords:** Software package, pipeline, intracranial recordings, group-level analysis, HGA, time-frequency decomposition, non-parametric statistics, Brainstorm.

**Graphical Abstract**

<figure><img src="image_2.png" /></figure>

# MIA: An Open-source Toolbox for Multi-patient Intracerebral EEG Analysis

## Pipelining the analysis of iEEG clinical signals to address neuroscientific questions

Electrode shaft implantations vary in anatomical location across patients

<figure><img src="image_3.png" /></figure>

Common sampled region
Here: left occipital gyrus

Functional activity identified across patients

<figure><img src="image_4.png" /></figure>

Group level activity computed, with corrected significance and consistency estimates

<figure><img src="image_5.png" /></figure>

Further details and funding information at www.neurotrack.fr/mia

Anne-Sophie Dubarry et al. 2021 @DubarryAnne

DOI to be completed
# 1. Introduction

Intracranial EEG (iEEG) is a method used for the investigation and treatment of various pathological conditions including drug-resistant epilepsy (Bancaud and Talairach, 1965). Used for clinical diagnosis, this inpatient procedure involves surgically implanting various electrodes, either on the neocortex (electrocorticography, ECoG) or intracerebrally (stereoelectroencephalography; SEEG), (Chauvel et al., 2019; Isnard et al., 2018; Minotti et al., 2018), then continuously recording the neural activity for about one to two weeks in an epilepsy monitoring unit (EMU). The recorded signals are examined to precisely locate and delineate the minimal pathological brain tissue or networks that need to be resected to render the patient seizure-free. At the same time, the procedure seeks to identify regions that need to be preserved in order to avoid post-surgical functional deficits (e.g. in language or memory abilities; Byrne, 2016; Lachaux et al., 2012; Trébuchon et al., 2020).

iEEG signals are recorded *in vivo* directly from cortical structures. They benefit from limited source mixing, including from artifactual sources such as eye blinks, and they arguably have very high signal-to-noise ratios (e.g., Ball et al., 2009). The information present in the signal can be revealed through various types of analyses, such as amplitude variations in time-locked averages (evoked responses in local field potentials), or the evolution of spectral content over time (time-frequency activities) (Tallon-Baudry and Bertrand, 1999). These analyses are used to provide fine-grained mappings of the epileptogenic zone and of cognitive networks to inform neurosurgical procedure decisions (e.g., Trébuchon & Chauvel 2016; Trébuchon et al., 2020).

In addition to the clinical assessment, intracranial recordings provide unique opportunities to perform cognitive neuroscience research, since patients can generally participate without interfering with their clinical procedures, under research protocols that are fully reviewed and approved by Institutional Review Boards. While the circumstances and population from which the signals are recorded are very specific, a variety of studies have linked the local and inter-regional activities to general properties of the cognitive processes involved in perception, attention, language, or memory (Flinker et al., 2018; Jacobs and Kahana, 2010; Lachaux et al., 2012; Llorens et al., 2011).

Patients with epilepsy may have had some cortical reorganization to compensate for epileptic dysfunctional tissue; therefore, for the signal analysis research to conclude on cognitive processes more broadly (i.e., on a population not suffering from epilepsy), two strategies are often used in such
analyses: (1) remove the signals recorded in pathological tissue, and (2) combine the results from multiple patients with varying epileptic sources (Flinker et al., 2018). With this second strategy, one major difficulty is the variability in the location of iEEG electrodes, which are implanted solely on clinical indications. While pharmaco-resistant epilepsy tends to concern certain brain regions and networks (e.g., temporal lobe epilepsy), the specific locations of electrode implantation may vary substantially from patient to patient. In turn, this anatomical variability may or may not be associated with variability in the recorded functional responses. Within a labelled anatomical region, various factors can contribute to the variability of the recorded responses: the inevitable approximations made during atlas co-registration and electrode labelling, actual heterogeneous functions within this region, or the potentially diverse influences of neighboring sources whose actual distance varies across patients.

These characteristics stand as a major challenge for performing group analysis and population wide inferences. They prevent the standard approaches used in non-invasive neurophysiology where a systematic spatial sampling is performed. In such a case, aggregating the data from multiple participants is achieved with various approaches such as standardizing the recording sensors placement (as in surface electroencephalography), or localizing the sources to alleviate cross-recordings variability of head position (as in magnetoencephalography).

Analysis tools dedicated to iEEG data were first developed to support clinical explorations. Some hospitals have developed routine pipelines which allow the analysis of task induced responses within a few hours (Cheung and Chang, 2012; Lachaux et al., 2007; Miller et al., 2007; Schalk et al., 2008). These so-called “localizers” or functional mapping procedures are developed in a concerted effort between the clinical staff and researchers/engineers to reach the clinicians' needs. The intended target is each individual patient (i.e., single patient analysis).

In cognitive neuroscience, the major software suites have been initially developed for non-invasive data (MEG or surface EEG): Brainstorm (Tadel et al., 2011), Fieldtrip (Oostenveld et al., 2011), MNE-Python (Gramfort et al., 2013), EEGLAB (Delorme and Makeig, 2004). More recently, some efforts have been deployed to make these tools compatible with iEEG data, notably the registration and labelling of contacts within the anatomy, or the automatic computation of various offline re-referencing montages (e.g., average, bipolar, local). However, no explicit procedure has been developed that would consider the difficulties induced by patient idiosyncrasies for group analysis.
Accordingly, it is not uncommon for iEEG empirical studies to implement the analysis of multiple patients using idiosyncratic methods (e.g. Anders et al., 2019; Barbeau et al., 2008; Dubarry et al., 2017; Nakai et al., 2017; Riès et al., 2017; Sahin et al., 2009), but this is done without introducing a systematic approach. To our knowledge, only two studies specifically addressed this issue: (Kadipasaoglu et al., 2014) presented a methodology based on surface-based co-registration and the use of mixed-effects multilevel analysis. More recently, (Magnotti, 2020) developed an analysis and visualization tool using R code that addresses multi patient iEEG data integration and which is fully standalone.

In this context, our goal was to design and implement an efficient open-source software for the analysis of multi-patient iEEG in the temporal and time-frequency domains. We developed an analysis pipeline that aims at capturing the steps and intuitions that an expert neurophysiologist would follow when combining single case patterns into a qualitative group description. At the same time, the approach had to go beyond intuitions to allow systematic and reproducible analysis. The MIA tool was developed and tested in the course of conducting various research projects (Anders et al., 2019; Dubarry et al., 2017; Trébuchon et al., 2020). We illustrate the use of the resulting toolbox using an SEEG dataset from a group of patients performing a language task (previously reported in Dubarry et al., 2017).
## 2. Toolbox architecture

MIA is an open-source toolbox developed under MATLAB (versions 2012a to 2017a; The Mathworks, Natick, MA), and therefore MIA can be installed on all hardware platforms and operating systems that MATLAB supports. Freely available for download under the GNU General Public License version 3 and accessible from our public Git repository (https://github.com/MIA-iEEG/mia), MIA is also available as a plugin for Brainstorm (release September 2021; Tadel et al. 2011). This open distribution and version control provides fully transparent content to motivate future collaborative work, such as developments of specific aspects of the tool to suit specific needs (see "Future developments" for examples). The general workflow within MIA consists in four five stages as illustrated on Figure 1 and detailed below.

<figure><img src="image_6.png" /><figcaption>Figure 1. MIA user general workflow and data processing stages.</figcaption></figure>

Following a standard clinical terminology (which may vary among sites, cf. Holdgraf et al., 2018), we will use the term *electrode* to refer to each of the brain probes, strips, or grids implanted in the patient's brain. In turn, each electrode typically includes between 5 and 15 distinct recording *contacts*
on a probe or strip, or up to 64 recording *contacts* (8 x 8) on a grid. The contacts are referenced using one of various *montages*, such as a *bipolar* montage (the difference between pairs of neighboring contacts) or a *monopolar* montage (the difference between each contact and a common reference). The signals recorded from these pairs of contacts are referred to as *channels*, and we note that channels of data can be arbitrarily re-montaged after the initial acquisition of the data.

## 2.1. Data preparation and import

The first MIA processing stage consists in importing the intracranial data (see Figure 1). MIA does not process raw signals. Therefore, the data require “pre-processing” steps which we describe here and that must be performed with external tools.

Data preparation can be achieved in Brainstorm (Tadel et al., 2011) where the MIA plugin is automatically downloaded from the menu “Plugins -> sEEG -> mia -> install”. MIA supports any of the input file formats supported by Brainstorm, including Deltamed, Micromed, Nihon Kohden, Nicolet, BrainAmp, ED, etc. The compatibility of MIA with these popular tools ensures that MIA can be used with virtually any file format currently available in clinical settings. It guarantees that the import process remains as robust as the process implemented in these well-established tools.

Optional pre-processing steps include:

* filtering (continuous data filtering is recommended),

* artifact rejection (marking segments on continuous data is recommended),

* bad channel rejection.

Irrespective of these or other steps, MIA being dedicated to the analysis of evoked neural activity, it is mandatory to perform an epoch segmentation around the stimulus of interest. These epochs should include enough temporal points to disregard edge effects (see section “Normalization and time-frequency extraction” for details). All these pre-processing steps should be documented as they are performed. It seems particularly relevant to highlight this point in the context of clinical data, which may be diverse across patients and sometimes collected from different research centers, etc. Performing the pre-processing steps should ensure that the dataset fed to MIA has a consistent structure across patients.
After data preparation, from Brainstorm, in processes menu ("Run" button), the custom process "File
-> Export data to MIA database" automatically imports the selected data into MIA and launches MIA
main GUI on the newly generated database.

As an alternative to Brainstorm, a legacy option is available if the data have been processed with BrainVision Analyzer (Brain Products GmbH, Munich, Germany). In this case the appropriate result to be used as input for MIA is composed of three files per patient (e.g., .eeg, .vmrk, .vhdr) containing all trials and which are read using external functions from EEGLAB (Delorme and Makeig, 2004). In this configuration, MIA can import multiple dataset at once, each containing the data from one patient and a default patient name is suggested for each. For scripting purposes (importing several patients at once), there is the option to exclude one or several contacts from each patient import by creating, in the patient directory, a bad.csv file containing the labels of the undesired channels (i.e., "bad channels"). The data from these channels will not be imported.

With the steps described above a MIA database is created which contains one folder per patient. Each of these folders will contain all the files that will be loaded at the later steps of the analyses. At the import stage (see Figure 1), one structured file is created (<patientID>_signal_LFP.mat; see Table 1) per patient, containing the three fields: "F", the imported signals in the form of a 3D matrix (Channels x Samples x Trials), "Labels": the channel labels and, "Time": a time vector.

In summary, these processing steps organize the sometimes-variable clinical data in a systematic data structure specifically designed for the follow-up analyses.

## 2.2. Data sanity check

After importing the data, MIA allows visualizing the imported data through an interface labeled "Sanity check" (see Figure 1). This graphical user interface displays general information about the recordings. The logic of this and other interface displays available throughout the data analysis pipeline is to keep the users close to their data, providing multiple opportunities for observing data features, and for performing visual controls at each of the processing steps. Of course, the data themselves (matrices) are also fully available if deeper controls are deemed necessary. The Sanity Check Display includes the number of trials, the total number of contacts (recording sites), the number of electrodes, and the number of contacts per electrode. The interface also shows the signals averaged over trials under two forms: a "butterfly" plot in which all time series are displayed along the same temporal axis, and a 2D
(raster plot) image on which the contacts' time series are stacked horizontally in alphabetical order and color codes for the amplitude.

At this stage, the user may visually identify and reject contacts that show artifactual activity on the averaged signals. Labeling a contact as “bad” at this stage will not permanently remove the contact but will store the information inside a channel description file (<patientID>_channels.tsv), following iEEG-BIDS recommendations (Holdgraf et al., 2018), which allows excluding the corresponding channel from subsequent analysis. Exclusions thus concern individual channels leaving all other signals unaffected. However, it is important to note that noise in one channel can affect all channels during the statistical procedure described later, where surrogate distributions are constructed by mixing channels. If such statistics are planned, then channel exclusion should be performed with said purpose in mind.

MIA allows the use of two montages: monopolar and bipolar. A monopolar montage consists in referencing all signals to a single recording site, typically with the hard-wired reference used during recording. This reference is either a scalp electrode (for instance Fz in a standard 10-20 surface montage), or a 'low variance' contact in iEEG (e.g., a contact located in the skull). A bipolar montage consists in subtracting the signals recorded from pairs of neighboring contacts, a procedure that is performed within each electrode. This montage allows capturing more focal activities by subtracting the (presumed distant) activity that is recorded equally in neighboring contacts. These two montages are available in MIA. It is important to note that the scheme used in the bipolar montage is exhaustive (e.g., c2-c1, c3-c2, c4-c3) rather than sparse (e.g., c2-c1, c4-c3). In the current version, if the data from a contact is missing, MIA will just follow the numbering order (e.g., if contact c2 is missing, MIA will create a bipolar channel c3-c1).
## 2.3. Normalization and time-frequency extraction

The user can choose to explore either the broadband signals or a specific frequency band (e.g., high gamma activity, HGA, between 80 and 150 Hz) (see Figure 1). In both cases, the data are normalized by a z-score procedure and the result is a 3D matrix of z-scored values.

If the broadband option is chosen, the z-score normalization procedure consists in subtracting at each time point the average over the baseline period and dividing by the standard deviation over this same period.

The time-frequency (TF) decomposition option is more elaborate and requires a detailed description, illustrated in Figure 2.

<figure><img src="image_7.png" /><figcaption>Figure 2. Extraction of time-frequency time courses. Input and output data are both 3D matrices (N Channels x N TimeSamples x N Trials). The algorithm presented here follows three steps: 1) Time-Frequency decomposition among several user-defined frequency bands; 2) z-score normalization against baseline to compensate for the 1/f distribution of power across the frequency spectrum; and 3) Averaging over normalized time-frequency band activities. The example shows frequencies chosen</figcaption></figure>
on a very large bandwidth to highlight the 1/f correction (see section “Time-frequency extraction” for further details).

In such decomposition, signal power modulations across time, frequency and trials are quantified by means of a continuous wavelet transform (CWT, Mallat, 2009), with Morlet wavelets as basis functions. This decomposition consists in convolving the neurophysiological signal s(t) with a Gaussian-shaped complex exponential function as described in (Bénar et al., 2009). The wavelet oscillates around the center frequency $f_0$, with a number of cycles ($\xi$) of the order of the oscillation parameter:

$$ \xi = 2\pi f_0 \sigma \qquad (1) $$

The length of the wavelet (d) varies accordingly. At 95% of the area under the Gaussian envelope the duration of the wavelet is given by

$$ d = 2 \frac{\xi}{\pi f} \qquad (2) $$

For instance, a 7-cycle wavelet transform corresponds to a wavelet duration of 55.7ms at 80Hz. The user defines the number of cycles $\xi$ (default = 7). The optimal number of cycles depends on the activity of interest: generally lower for low-frequency activity, and higher for high frequency oscillations. Under this approach, spectral bandwidth increases with frequency.

To compensate for the 1/f distribution of power across the spectrum, i.e. decreasing power with increasing frequency (Buzsáki, 2006; Milstein et al., 2009), the decomposition is first performed at several frequencies within the range of interest (e.g. 80-150 Hz). The space between the frequency bins within the range of interest is defined by the parameter “step” (e.g., 80, 90, 100, ..., 150 Hz for a 10 Hz step between 80 and 150 Hz). At each frequency ($f_i$), the magnitudes of the wavelet coefficients ($|w|$; module of the complex coefficients) are z-scored against their respective baselines ($b$). Using the pre-stimulus period as baseline period (rather than the entire epoch) ensures that the activity of interest is not reduced (Roehri et al., 2016). Finally, all z-scored values are averaged across the $n$ frequencies within the band:

$$ wn(t) = \frac{1}{n} \sum_{f} \frac{|w_{f_i}(t)| - |w_{f_i}(b)|}{|\sigma_{f_i}(b)|} \qquad (3) $$
For examples of a similar approach, although implemented with a Hilbert transform, see (Hamamé et al., 2014; Vidal et al., 2012, 2010) among many others. We apply this procedure (either normalization only or time-frequency and normalization) on each channel and each trial to obtain N<sub>Channels</sub> x N<sub>Tim</sub> eSamples x N<sub>Trials</sub> matrix ("Output" in Figure 2).

Using this approach, activity that is phase locked to the stimulus across trials (i.e., evoked) and activity that has temporal jitter across trials (i.e., induced) will result in similar increases of energy in the time-frequency domain (Tallon-Baudry and Bertrand 1999). An option to subtract the evoked component before computing the induced component is available.

## 2.4. Statistical analysis

For each patient, we assess the statistical significance of post-stimulus signal variation versus the baseline, by performing a two-stage statistical procedure: 1) a t-test and 2) a correction for multiple comparisons.

First, a one-sample Student t-test against 0 is computed on the z-scored data across trials at each channel and each time point of the data. The result is a t-value statistical map (N<sub>Channels</sub> x N<sub>TimeSamples</sub>), which is converted into a probability map (p-map) based on the known distribution of the t-statistic. This conversion depends on the number of degrees of freedom (i.e., number of trials minus 1). This probability activation map is then thresholded to a pre-established alpha value (default at p = 0.001; see example in Figure 3). This thresholded map is *not yet* corrected for multiple comparisons and should be interpreted cautiously; however, it is displayed to help appreciate the impact of the multiple comparison step described in the next paragraph.

With the above procedure, the t-test is repeated as many times as there are time samples and contacts, resulting in the well-known multiple comparison problem. To solve this problem in the time-domain, MIA combines two established procedures, i) a minimum duration of significance threshold (similarly to temporal clusters in Maris and Oostenveld, 2007) and (ii) a bootstrap procedure that enables testing only one condition versus zero (Delorme and Makeig, 2004). Thus, trials from the original data set are randomly selected (with repetitions allowed, Efron, 1992) and periods with significant activity within the baseline window are identified (i.e., p < a; default a = 0.001; Table 1). The procedure is repeated N times (default N = 500; Table 1), while the maximum number of contiguous points passing the significance threshold across all contacts at each iteration is pooled into a surrogate
distribution of maximum durations (see Figure 3). The 95% quantile of this bootstrap distribution provides the significance duration threshold. Finally, the statistical map is thresholded using the duration threshold obtained with the bootstrap procedure; in other words, clusters identified in the first step must be longer than the bootstrap-determined threshold This procedure (taking the maximum at each iteration) allows alleviating the multiple comparison problem jointly across time points and contacts.

As noted under “Data Sanity Check”, a channel with a large amount of noise during the baseline that is not excluded at the pre-processing stage could perturb the statistical analysis. At each bootstrap iteration, the artifacted channel will show the cluster with maximum duration of activation that will be kept for the surrogate distribution. This will shift the distribution to unduly high values. The duration threshold would then be over-estimated for the whole set of channels, resulting in an over-estimation of the threshold for detecting post baseline activations, many of which would go undetected. Once again, this would apply to all contacts and not only the noisy channel itself. Artifacted channels must either be rejected during the pre-processing or be marked as “bad channels” during the sanity check.

As a final result of this processing stage, MIA provides a user interface displaying three maps side by side: the map without threshold (z-scores), the map thresholded but uncorrected (uncorrected t-values), and the map thresholded and corrected (corrected t-values). These displays illustrate the different steps of the procedure and allow the comparison of their outputs. With this flexible graphic interface, the user gets a broad view of the data under study. In addition, all the parameters of the analysis (e.g., *alpha*, *p*, *iterations*; Table 1), the resulting data structure, and the display (selected channels) are saved for future reference and ease of reuse and reproduction.
<figure><img src="image_8.png" /><figcaption>Figure 3. Statistical analysis pipeline. On the pre-processed epoched data (3D matrix; "Input"), MIA computes one sample Student t-tests across trials at each channel and each time point of the data ("Operation"). The result of this operation is an uncorrected statistic map ("Example"). To address the multiple comparisons problem across timepoints and contacts, MIA estimates a minimum duration threshold for consecutive t-values with a non-parametric approach. Finally ("Output"), periods of significance which do not exceed this threshold are removed from the statistical map.</figcaption></figure>

## 2.5. Group analysis

As a final processing stage, MIA allows performing group analyses (see Figure 1). A "study" can be created from multiple patient data that have similar analysis parameters: frequency bands, montage, etc. The analysis consists in grouping each patient's contacts by anatomical region, thus effectively switching from a database of patients to a database of regions.

* **Labeling table**

This grouping procedure requires a labeling table which contains, for each patient, the correspondence between the iEEG recording channels and anatomical regions. This table might be obtained from the
clinical procedure, for example from the observation of fused CT-MRI images by an expert. Alternatively, it might be constructed with automated tools (Blenkmann et al., 2017; Deman et al., 2018; Groppe et al., 2017; Medina Villalon et al., 2018), which integrate the whole pipeline of CT-MRI co-registration and contact labeling. Currently, MIA supports labeling table constructed with Brainstorm, Fieldtrip or a custom table organized such that each line describes one recording contact, with four simple fields across four columns: patient ID, contact label, laterality (L or R) and region anatomical label. The region labels can be arbitrarily decided by the user, but they must be consistent across patients.

## Inclusive/exclusive

For monopolar montages, the correspondence between contacts and regions is provided directly in the labeling table. For bipolar montages, a decision needs to be made when the two (adjacent) contacts of a bipolar channel are classified in two different regions in the labeling table. The user can choose to label a channel as belonging to a region either: (A) if both of the contacts do belong to the region ("exclusive" rule) or (B) if at least one of the two contacts belongs to the region ("inclusive" rule). By default, bipolar channels for which at least one of their contacts was labeled as being outside the brain, or when both contacts are located in white matter, are rejected from subsequent analysis.

<figure><img src="image_9.png" /><figcaption>Figure 4. Intracerebral stereotaxic EEG (SEEG) recording contacts and regions, visually labeled based on the atlas (Mai et al., 2008). Left panel: a slice of the Mai et al atlas (permission pending). The various regions are represented with different colors. Middle panel: for one patient, the corresponding MRI slice in which an electrode is visible; the contact labels are written near the regions (Fop: Frontal</figcaption></figure>
operculum; WM: White matter; PrG: Pre-central Gyrus; out: outside of cortex). Right panel: an example section of a labeling table containing four columns: patient, contact, laterality and region (PT_10 is the patient displayed on the middle panel).

* **Correlations (channels/patients)**

To assess the consistency of the activity recorded within a region, MIA computes two correlation measures, the inter-channel correlation (across the time courses of all relevant channels within the region), and the inter-patient correlation (across the time-courses averaged per patient within the region). Both correlation metrics are defined as the mean of all pair-wise Pearson correlations within a region, either among contacts or among patient average time-courses. The MIA interface does not provide a test of statistical significance for these correlations. An additional standalone script is provided to perform a bootstrap test which is implemented as follows (*mia_get_roi_permute.m*). For every region of interest, the structure of patient and channel sampling available to compute the observed correlation is preserved to draw random samples of data from among all available regions. A correlation is computed for each of these draws, and values are pooled into a surrogate distribution of correlations (corresponding to the null hypothesis of no relationship). A threshold value is found from the surrogate distribution (cumulative histogram above 1-p), in order to determine if the naturally observed correlation is deemed significant or not.

* **Sign flip**

In the particular case of evoked local field potentials computed in the time domain using a monopolar montage, differences in sign may result from differences in the relative positions of the generator and the recording contact (Buzsáki et al., 2012; Lopes da Silva and Van Rotterdam, 2005). This feature is often referred to as polarity reversal. It can be highly informative when examined closely but, in a group analysis, it can potentially compromise the interpretability of any averaging performed across contacts or patients. To circumvent it, we implemented a procedure intended to capture a neurophysiologist's interpretation of signed potentials and their aggregation.

The procedure is based on the idea of compensating for polarity reversals by flipping the sign of some of the signals before averaging across contacts and patients. The challenge is in deciding whether and which of the signals should be flipped. Our rationale was to identify the time course that is the most negatively correlated to all the other time-courses. This time course (and those that are most positively correlated with it) are then flipped to better match the polarity of the remaining channels.
More concretely, within each region, the Pearson's correlation matrix is computed between all channels. This correlation matrix is thresholded for negative values (default threshold value -0.7). In general, high negative correlations that lie beyond the threshold value are considered to result from polarity reversals and to require sign flipping. The number of negative correlation values passing the threshold is counted and the channel with the highest number is identified as an anchor for the polarity reversal. All channels positively correlated (above threshold) with this "anchor" have their sign flipped. Conversely, correlations that lie below the threshold parameter (e.g. between -0.7 and 0.7) are considered insufficient to reflect a polarity reversal that would require sign flipping. This approach was preferred over the computation of absolute values to preserve the shape of the signal and the information carried by inversions of electric potential.

* Show only significant

There is the option to limit the analysis and displays to the contacts for which the post-stimulus signal variation versus the baseline is significant (corrected for multiple comparison, see above). If this option is used, the analysis is performed only with contacts that show at least one period of significant signal in the time window of interest.

## 2.6. Additional features

* Grouped raster plots

To explore the single trial data within a group analysis, we developed a raster plot viewer which allows the user to display, for each region, the single-trial data in a stacked layout (Delorme et al., 2015). The stack is organized by electrodes and patients. When needed, trials can be ordered by a parameter chosen by the user (e.g., response times, trial order, etc.)

* Export results by region(s)

From the "Group analysis" GUI the user can export one or several region(s) to the Matlab workspace. This gives the user the opportunity to analyze the data from regions with his/her own custom code/software. This export creates a variable inside the Matlab workspace called "roi" (region of interest; a cell array which dimension is 1 x number of selected regions) containing various variables concerning the region(s) (e.g., timeseries, labels, correlations)
<table><caption>Table 1: List of processing parameters ordered by processing stage.</caption><thead><tr><th>Processing stage</th><th>Parameter</th><th>Description</th><th>Default values</th></tr></thead><tbody><tr><td rowspan="3">Import</td><td>OPTIONS.outdir</td><td>Output directory</td><td>Prompt user at first execution</td></tr><tr><td>OPTIONS.subjects</td><td>List of patients ID</td><td>"</td></tr><tr><td>sInput</td><td>Input file names</td><td>"</td></tr><tr><td rowspan="10">Time-Frequency decomposition</td><td>Output filename</td><td>&lt;patient&gt;_signal_LFP.mat</td><td></td></tr><tr><td>OPTIONS.outdir</td><td>Output directory</td><td>MIA database path</td></tr><tr><td>OPTIONS.subjects</td><td>List of patients to process</td><td>"</td></tr><tr><td>OPTIONS.mtg</td><td>Montage</td><td>'Monopolar'</td></tr><tr><td>OPTIONS.modetf</td><td>Time-Frequency method used to decompose signal</td><td>'Morlet'</td></tr><tr><td>OPTIONS.freqs</td><td>Decomposition frequencies</td><td>[0:0:0]</td></tr><tr><td>OPTIONS.removeEvoked</td><td>Option to remove evoked response (0 or 1)</td><td>0</td></tr><tr><td>OPTIONS.zbaseline</td><td>Baseline period used for normalization</td><td>Available time window without edge effects</td></tr><tr><td>OPTIONS.overwrite</td><td>Overwrite existing results ('Yes' or 'No')</td><td>'yes'</td></tr><tr><td>Output filename</td><td>&lt;patient&gt;_&lt;montage&gt;_&lt;method&gt;_data_LFP.mat</td><td></td></tr><tr><td rowspan="5">Statistics</td><td>OPTIONS.outdir</td><td>Output directory</td><td></td></tr><tr><td>OPTIONS.subjects</td><td>List of patients to process</td><td>"</td></tr><tr><td>OPTIONS.alpha</td><td>Alpha value</td><td>0.001</td></tr><tr><td>OPTIONS.nboot</td><td>Number of bootstraps</td><td>500</td></tr><tr><td>OPTIONS.baseline</td><td>Baseline period used for statistics</td><td>[-0.8, -0.05];</td></tr><tr><td>Output filename</td><td>&lt;patient&gt;_&lt;montage&gt;_&lt;method&gt;_stats_LFP.mat</td><td></td></tr><tr><td rowspan="4">Grouping</td><td>OPTIONS.twin</td><td>Time window of interest</td><td>[-0.5, 1.5]</td></tr><tr><td>OPTIONS.loctable_name</td><td>Labeling table file path</td><td>"</td></tr><tr><td>OPTIONS.inclusive</td><td>Inclusive/Exclusive rule for labeling bipolar channel</td><td>0</td></tr><tr><td>Output filenames</td><td>m_table_&lt;Table Name&gt;.mat<br>&lt;study&gt;_&lt;Npatients&gt;_&lt;method&gt;_&lt;montage&gt;.mat</td><td></td></tr></tbody></table>

## 3. Results

The processing options in MIA are illustrated using the data reported in (Dubarry et al., 2017). The original goal of the study was to assess the degree to which word retrieval involves parallel (simultaneous) processing across brain regions and hence, presumably, cognitive processes. To address this question, 15 patients with epilepsy implanted with SEEG for pre-surgical diagnosis purposes were asked to perform a picture-naming task. Intracerebral EEG data were recorded during the task. The signal was time locked to stimulus onset. All data are publicly shared for transparency and reproducibility on the Zenodo open-access repository (DOI: 10.5281/zenodo.4767855)
### 3.1. Import data and sanity check

For the standalone version of MIA, the user has to add the package to Matlab path (by using the command `addpath`). Then entering the 'mia' command in Matlab prompt opens the main interface (Figure 5).

<figure><img src="image_10.png" /><figcaption>Figure 5. Main GUI window. (A) The database: a button allows browsing for a database directory. The full path of this directory is displayed. This directory is initially empty and will be filled by MIA as the user analyzes the data. (B) The left column shows the list of imported patients, (C) The middle column shows the processed data (either z-scored or TF decomposed), sorted by patient. The lower box shows the parameters of the statistics used to compute the selected data file. (D) The right column shows the group analysis data processed (upper box: data grouped over multiple patients; lower box: anatomical labeling table). At the first execution, the user will be prompted to set a database folder. This path can be changed later (e.g. for subsequent analyses of other datasets) by pressing the "Database" button.</figcaption></figure>

For the Brainstorm plugin MIA version, the user exports the data into a MIA database and the main interface automatically opens and sets the database folder to the newly generated one.
The left panel (Patients) contains a list of the patients' labels in the current database. Each entry of the list corresponds to a new patient dataset which has been imported (button "import"). After selecting a patient from the list, the user can display the sanity check GUI by pressing the "Display" button (Figure 5).

<figure><img src="image_11.png" /><figcaption>Figure 6. Sanity Check GUI window. (A) General information about the dataset: Patient identification, number of trials, total number of contacts, total number of electrodes and number of contacts per electrode. (B) Average time series shown on a butterfly view (amplitude of all contacts on a shared Y-axis). (C) Average time series shown on an image view (electrodes on the Y-axis, amplitude on color a scale)</figcaption></figure>

## 3.2. Time-frequency and statistical analysis

To normalize (z-score) the signals and/or extract time-frequency activity the user must press the "Frequency" button. A window will prompt the user to fill in the frequency range, its steps, the method of extraction (z-score normalization or Morlet Time-frequency decomposition), the number of cycles
(oscillation parameter; defined in equation 1) for Morlet wavelets, the baseline range, and the
montage.

To visualize the results of the time-frequency decomposition (or normalized broadband signals) and
the statistical computations, the user must select a data file and press the “Display” button located in
the Data panel. The resulting display is presented on Figure 7.

<figure><img src="image_12.png" /><figcaption>Figure 7. Single patient data and statistics display window. Results from one patient processed with the MIA pipeline. The three panels show channel stackplots. From left to right: 1) the averaged trials z-score (see Figure 3 output), 2) the result of the t-test against zero (uncorrected) and 3) the corrected t-values (bootstrap correction).</figcaption></figure>
Patient ID and data information are displayed at the top of the window. A dropdown menu shows the different parameters that were used to compute statistics on the datafile in view. If no statistics have been computed on the current dataset, the user will receive a warning pop-up message before the window opens.

Channels can be selected from the "Select Channel" button. This selection applies online for the three panels at the same time. It also applies for the plot of the z-scored time courses which is accessible from the button located underneath the left (z-scores) panel ("Plot Time Series" button). The selection of a set of channels sticks at the re-execution of the display: the same subset of channels will be displayed if this dataset is displayed again, even if the program has been closed in the meantime. Amplitude scales can be adjusted manually, separately for the z-scored values (lower left corner) and for the t-values (lower right corner). By default, these scales are set to the maximal value of the data. At the top of the window, a menu bar contains six buttons described in Table 2.

<table><caption>Table 2: Buttons from statistics display</caption><thead><tr><th>Button</th><th>Function</th><th>Comment</th></tr></thead><tbody><tr><td>Zoom in and Zoom out</td><td>Give access to the Matlab zoom feature</td><td>Applicable on any of the three panels (the two other panels will accordingly adapt to the same zoom accordingly)</td></tr><tr><td>Data cursor</td><td>Get the precise cursor information (time and amplitude)</td><td>Click anywhere on one of the three panels</td></tr><tr><td>JPEG button</td><td>Saves a snapshot of the window</td><td>A default "JPEGs" folder is automatically created in the main database directory</td></tr><tr><td>Pan tool</td><td>Pan the axes views</td><td>A native Matlab feature</td></tr><tr><td>DISP button</td><td>Hide or show the left and middle panels</td><td>These correspond to z-scores and uncorrected t-values</td></tr></tbody></table>

### 3.3. Group analysis

After creating a study and loading a labeling table, the user must press the display button from the Studies panel (Figure 5D) to open the group analysis GUI.
<figure><img src="image_13.png" /><figcaption>Figure 8. Group analysis GUI window. Interactive Regions Of Interest (ROI) table for group analysis. (A) At the top of the GUI is indicated the name of the labelling table which was used for the group study. (B) The user can switch between studies. (C) The interactive table shows the list of ROIs for the current study; the content of the table can be sorted according to any of the fields by clicking on the column header. (D) Various types of displays can be used for activities averaged by region(s). (E) Single trials activity can be displayed in a stackplot view ("Rasters"); if response times (RT) are available, trials are automatically organized by patient, channels and RT (see Figure 10). (F) There is the option to use sign flip and show only significant (see section on Implementation) and to export one or several selected region(s).</figcaption></figure>

This GUI contains an interactive table of regions where the data can be sorted by different fields (number of patients, contacts, anatomical regions, channel correlation, patient correlation) and provides different display options. From the interactive table, the user can select one or more regions and apply various display modes. The button "Show table" displays the interactive table. The user can select one or several regions at a time and apply two types of displays: averaged activities (Display ROIs panel) and single trial activities (Rotor Plot panel).
For averaged activities, three main types of display are available (Figure 8, under the box “display ROIs”). The first type, “ROIs”, displays the time series of the contacts which belong to the region(s) selected in the bottom list. The top panel of the ROIS display shows signals averaged per patient, while the lower panel shows the mask of significance for all contacts, stacked on top of one another. The second type, “Mean ROI”, displays the time series averaged across contacts and patients for the selected region(s) in the three possible representations, namely a “stacked plot” (Figure 9A), a “butterfly plot” (Figure 9B) or a “flat display” (Figure 9C) in which the averaged activities are displayed as an image and the color codes for the amplitude (z-score); the regions are stacked one over the others along the x-axis. In the two latter displays, stacked and flat, the signals are ordered vertically along the regions' first onset latency (i.e., the timing of the first sample exceeding 1.96 z-score). Finally, the third type of display that can be selected, “One patient data” shows the patient's average per region, for one or several patients.

<figure><img src="image_14.png" /><figcaption>Figure 9. Multi-patient averages for different ROIS presented under three different displaying options: (A) “Stacked plot” display (B) “Butterfly plot” (C) “Flat display”</figcaption></figure>

For single trial activities, the user has two main display options: 1) a vertically stacked single trials image view below the averaged activities, shown with the “ROIs” button (Figure 8). 2) a display of the vertically stacked single trials image only (recommended for regions involving many contacts).
In order to display single trial activities, the user has to load the signals by using the "Load Rasters" button. This step was implemented to reduce the time required to open the group interface (as it may take time to load all of the single trial signals). The button "Get RT file" allows the user to load a file that must contain at least the following columns: "Patient", "Order", "RT", "InData". The patient and RT column are used to display the single trials ordered by response time and to indicate the RT by a black point. The Order and InData are intended to capture, respectively, a trial order (in case the RT is not used) and the trials which would be present in the behavioral data file ("RT file"), but not in the data (e.g. because they were rejected from the data during the pre-processing).
<figure><img src="image_15.png" /><figcaption>Figure 10. Single trials raster plot for region "L.TTG1" (left anterior transverse temporal gyrus). The data from five bipolar channels located in the primary auditory cortex in two patients (PT_16 and PT_22). Trials are ordered by response times, activity is locked to the response. The increase following the response is presumably related to hearing own speech.</figcaption></figure>
# 4. Discussion

We introduce MIA, a new toolbox for performing group studies of multi-patient iEEG data in the time and frequency domains. While our main objective is to facilitate the analysis of intracranial recordings in the context of cognitive neuroscience research, it is also possible to use it for exploring single patient data in more clinical settings. The pipeline captures the most common steps that a neurophysiologist might follow to assess and group the data of various patients. For group analysis, MIA allows the user to handle inter-patient variability in electrode position and time course of both broadband and time-frequency activities. It provides comprehensive summaries of the spatio-temporal dynamics of brain activity observed across patients, with statistical tests both within condition (i.e., against a pre-stimulus baseline) and across conditions (e.g. two different levels of a cognitive factor).

The most central and original aspects of the proposed analysis pipeline are (1) the normalization of the signal power at each frequency band, which substantially reduces the 1/f effect (for a thorough discussion on the issue see Haller et al 2018), (2) the consideration given to the multiple comparison problem (MCP) along the temporal dimension and across all contacts from a single individual (based on the proposals by (based on the proposals by Efron, 1992; Maris and Oostenveld, 2007), (3) the interactive tables and visualization interface, specifically adapted for handling multiple patient data at the level of grand averages, channel averages, and single trials, and (4) the open dissemination of a version code, which contributes to the transparency and reproducibility of the analysis, even for those who are not familiar with handling large amounts of source code.

We highlight other features that may require explicit decisions from users deploying their own analyses. One of the implemented analyses allows determining whether post-stimulus signal activity differs significantly or not from the corresponding baseline, thus yielding the channel labels "significant" or "not significant". First, in the non-parametric approach, which account for the MCP, we advocate that both the baseline and the time window of interest have the same length so that the probability to obtain a cluster by chance is the same in the two time windows (on this issue see Piai et al., 2015b, 2015a). Second, in subsequent analysis (e.g. averages or correlations across contacts or patients), there is an option to include all contacts at stake (e.g. within a given region) or only those contacts that show significant activity within that region. Selecting a data subset on the basis of one test before performing another test on the same data is understood as "double-dipping", a practice that may distort the results to variable extents (the argument has been particularly developed for functional brain imaging: Button, 2019; Kriegeskorte et al., 2009). Users should be aware of links
between their selection criteria and the subsequent statistical tests to assess the interpretability of
their findings (for a thorough treatment, see supplementary materials of Kriegeskorte et al., 2009).

The time-frequency analysis has been programmed with a lot of flexibility in the choice of parameters
such as frequency bands, steps, time-windows etc. The choice of these parameters must be guided by
the research questions. For example, collapsing frequencies between 10 and 120 Hz may have little
interpretative value, as that may collapse different neural processes (Buzsáki, 2006). As another
example, choosing narrow time-windows may compromise the analysis of lower frequency bands
because of border effects. In addition, we highlight the use of the labeling table to code for the
anatomical classification of contacts; however, users are free to adopt other classification criteria (e.g.
functional responses in another data set) to group the contacts in the group analysis. The current
initiatives in developing large databases (e.g. OpenNeuro at openneuro.org) opens up perspectives for
processing large amounts of data and call for the possibility to test other types of classifications.

MIA was designed and developed with the idea of improving current practice in the group analysis of iEEG data. The program can input a wide variety of data formats through its connection with two well-established tools: Brainstorm and BrainVision. It includes rigorous descriptive and inferential signal processing and statistical methods that are explicitly described. Furthermore, the version-controlled open source format allows these methods to be edited, adapted or upgraded as needed in future developments, for example, alternative time-frequency decomposition methods such as Hilbert transform (Hamamé et al., 2014; Vidal et al., 2012, 2010) or Morse wavelet (Gadhoumi et al., 2012; Lilly and Olhede, 2010), or alternative re-referencing contact montages such as in (Mercier et al., 2017). Finally, MIA provides the user with as much information as possible regarding the data, at every step of the analysis, to discourage a potential “black-box approach” to data analysis, and the detailed steps of the analysis can be explicitly saved and/or scripted for full reliability and reproducibility.
## 5. Future developments

To perform multi-patient intracranial analysis within an atlas of interest, MIA requires an anatomical labeling table. Various tools have been designed for the automatic labeling of contacts (e.g. Groppe et al., 2017) within existing atlases (e.g "MarsAtlas" Auzias et al., 2016). In the current state, MIA supports labeling tables generated by Brainstorm, Fieldtrip, and custom tables created manually (e.g., from the expert visual inspection of anatomical locations). Additional reading functions will be developed promptly to ensure compatibility with other software (e.g. Gardel, Medina Villalon et al., 2018). In addition, the recent development of automatic labeling of all SEEG/ECOG contacts within any of the available atlases in Brainstorm opens the perspective to transfer quantitative information computed in MIA back into Brainstorm, for example correlation values, numbers of patients or numbers of contacts per region. The user would thus benefit from the ergonomic Brainstorm graphical tools to explore MIA results.

MIA is currently focused on measures of event related local activity and their variability across patients. Similar estimates of variability across different networks, involving more than one region and measured from event-related or spontaneous activities, will allow investigating other aspects of cognitive dynamics (Fries, 2005). Also, in the pipeline itself, we plan to develop modules for the processing of single trial features, such as the estimated onset time of a contact or region's significant activity (Coon and Schalk, 2016).

Finally, MIA supports both SEEG and ECOG data, although it was not tested with ECOG data; we plan to thoroughly test the toolbox with an ECOG dataset. And, because the challenges of interindividual variability are not restricted to iEEG, we also plan to combine the iEEG time series with source reconstructed data from MEEG signals. This would be particularly useful for simultaneous recordings of MEG/EEG and iEEG (Dubarry et al 2014, Badier et al 2019).

## 6. Conclusion

Cognitive neuroscientists increasingly rely on iEEG data to address scientific questions. The many sources of variability that may affect iEEG data pose a serious and largely unresolved methodological challenge for group analyses. To date, no consensual analyses procedure has emerged. We present a user-friendly tool, developed in collaboration with iEEG neurophysiologists and clinicians, that implements their basic steps for single case analysis and procedures for more elaborate multi-patient integration, with multiple parameter settings, and an open-source code.
Acknowledgements

We would like to thank Anais Llorens for sharing the dataset we used to illustrate MIA's functions; Vitoria Pia, Chotiga Pattamadilok, Manuel Mercier, and John Mosher for providing valuable feedback on previous versions of the manuscript; Joey Hsu, Shuai Wang, and Matteo Mascelloni for their contributions to the MIA online repository; and Francois Tadel for his guidance and recommendations regarding the toolbox.

References

Anders, R., Llorens, A., Dubarry, A.-S., Trébuchon, A., Liégeois-Chauvel, C., Alario, F.-X., 2019. Cortical Dynamics of Semantic Priming and Interference during Word Production: An Intracerebral Study. Journal of Cognitive Neuroscience 31, 978–1001.
https://doi.org/10.1162/jocn_a_01406

Auzias, G., Coulon, O., Brovelli, A., 2016. *MarsAtlas*: A cortical parcellation atlas for functional mapping. MarsAtlas. Hum. Brain Mapp. 37, 1573–1592. https://doi.org/10.1002/hbm.23121

Barbeau, E.J., Taylor, M.J., Regis, J., Marquis, P., Chauvel, P., Liegeois-Chauvel, C., 2008. Spatio-temporal Dynamics of Face Recognition. Cerebral Cortex 18, 997–1009.
https://doi.org/10.1093/cercor/bhm140

Bénar, C.G., Papadopoulo, T., Torrésani, B., Clerc, M., 2009. Consensus Matching Pursuit for multi-trial EEG signals. Journal of Neuroscience Methods 180, 161–170.
https://doi.org/10.1016/j.jneumeth.2009.03.005

Blenkmann, A.O., Phillips, H.N., Princich, J.P., Rowe, J.B., Bekinschtein, T.A., Muravchik, C.H., Kochen, S., 2017. iElectrodes: A Comprehensive Open-Source Toolbox for Depth and Subdural Grid Electrode Localization. Front. Neuroinform. 11. https://doi.org/10.3389/fninf.2017.00014

Button, K.S., 2019. Double-dipping revisited. Nat Neurosci 22, 688–690.
https://doi.org/10.1038/s41593-019-0398-z

Buzsáki, G., 2006. Rhythms of the Brain. Oxford University Press.
https://doi.org/10.1093/acprof:oso/9780195301069.001.0001

Buzsáki, G., Anastassiou, C.A., Koch, C., 2012. The origin of extracellular fields and currents--EEG, ECoG, LFP and spikes. Nat Rev Neurosci 13, 407–420. https://doi.org/10.1038/nrn3241

Cheung, C., Chang, E.F., 2012. Real-time, time-frequency mapping of event-related cortical activation. J. Neural Eng. 9, 046018. https://doi.org/10.1088/1741-2560/9/4/046018

Coon, W.G., Schalk, G., 2016. A method to establish the spatiotemporal evolution of task-related cortical activity from electrocorticographic signals in single trials. Journal of Neuroscience
Methods 271, 76–85. https://doi.org/10.1016/j.jneumeth.2016.06.024

Delorme, A., Makeig, S., 2004. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. Journal of neuroscience methods 134, 9–21. https://doi.org/10.1016/j.jneumeth.2003.10.009

Delorme, A., Miyakoshi, M., Jung, T.-P., Makeig, S., 2015. Grand average ERP-image plotting and statistics: A method for comparing variability in event-related single-trial EEG activities across subjects and conditions. Journal of Neuroscience Methods 250, 3–6. https://doi.org/10.1016/j.jneumeth.2014.10.003

Deman, P., Bhattacharjee, M., Tadel, F., Job, A.-S., Rivière, D., Cointepas, Y., Kahane, P., David, O., 2018. IntrAnat Electrodes: A Free Database and Visualization Software for Intracranial Electroencephalographic Data Processed for Case and Group Studies. Front. Neuroinform. 12, 40. https://doi.org/10.3389/fninf.2018.00040

Dubarry, A.S., Llorens, A., Trébuchon, A., Carron, R., Liégeois-Chauvel, C., Bénar, C.G., Alario, F.X., 2017. Estimating Parallel Processing in a Language Task Using Single-Trial Intracerebral Electroencephalography. Psychological Science 28, 414–426. https://doi.org/10.1177/0956797616681296

Efron, B., 1992. Bootstrap Methods: Another Look at the Jackknife, in: Kotz, S., Johnson, N.L. (Eds.), Breakthroughs in Statistics: Methodology and Distribution. Springer New York, New York, NY, pp. 569–593. https://doi.org/10.1007/978-1-4612-4380-9_41

Flinker, A., Piai, V., Knight, R.T., 2018. Chapter 43. Intracranial electrophysiology in language research, in: The Oxford Handbook of Psycholinguistics. p. 17.

Gadhoumi, K., Lina, J.-M., Gotman, J., 2012. Discriminating preictal and interictal states in patients with temporal lobe epilepsy using wavelet analysis of intracerebral EEG. Clinical Neurophysiology 123, 1906–1916. https://doi.org/10.1016/j.clinph.2012.03.001

Gramfort, A., Luessi, M., Larson, E., Engemann, D.A., Strohmeier, D., Brodbeck, C., Goj, R., Jas, M., Brooks, T., Parkkonen, L., H??m??l??inen, M., 2013. MEG and EEG data analysis with MNE-Python. Frontiers in Neuroscience 7, 1–13. https://doi.org/10.3389/fnins.2013.00267

Groppe, D.M., Bickel, S., Dykstra, A.R., Wang, X., Mégevand, P., Mercier, M.R., Lado, F.A., Mehta, A.D., Honey, C.J., 2017. iELVis: An open source MATLAB toolbox for localizing and visualizing human intracranial electrode data. Journal of Neuroscience Methods 281, 40–48. https://doi.org/10.1016/j.jneumeth.2017.01.022

Hamamé, C.M., Alario, F.-X., Llorens, A., Liégeois-Chauvel, C., Trébuchon-Da Fonseca, A., 2014. High frequency gamma activity in the left hippocampus predicts visual object naming performance. Brain and language 135, 104–14. https://doi.org/10.1016/j.bandl.2014.05.007

Holdgraf, C., Appelhoff, S., Bickel, S., Bouchard, K., D'Ambrosio, S., David, O., Devinsky, O., Dichter, B., Flinker, A., Foster, B., Gorgolewski, C., Groen, I.I.A., Groppe, D., Gunduz, A.,
Hamilton, L., Honey, C.J., Jas, M., Knight, R., Lachaux, J.-P., Lau, J.C., Lundstrom, B.N., Lee-Messer, C., Miller, K., Ojemann, J.G., Oostenveld, R., Piantoni, G., Petridou, N., Pigorini, A., Pouratian, N., Ramsey, N., Stolk, A., Swann, N.C., Tadel, F., Voytek, B., Wandell, B., Winawer, J., Zehl, L., Hermes, D., 2018. BIDS-iEEG: an extension to the brain imaging data structure (BIDS) specification for human intracranial electrophysiology.
https://doi.org/10.31234/osf.io/r7vc2

Jacobs, J., Kahana, M.J., 2010. Direct brain recordings fuel advances in cognitive electrophysiology.
Trends Cogn Sci 14, 162–171. https://doi.org/10.1016/j.tics.2010.01.005

Kadipasaoglu, C.M., Baboyan, V.G., Conner, C.R., Chen, G., Saad, Z.S., Tandon, N., 2014. Surface-based mixed effects multilevel analysis of grouped human electrocorticography. NeuroImage 101, 215–224. https://doi.org/10.1016/j.neuroimage.2014.07.006

Kriegeskorte, N., Simmons, W.K., Bellgowan, P.S.F., Baker, C.I., 2009. Circular analysis in systems neuroscience: the dangers of double dipping. Nat Neurosci 12, 535–540.
https://doi.org/10.1038/nn.2303

Lachaux, J.-P., Axmacher, N., Mormann, F., Halgren, E., Crone, N.E., 2012. High-frequency neural activity and human cognition: past, present and possible future of intracranial EEG research.
Progress in neurobiology 98, 279–301. https://doi.org/10.1016/j.pneurobio.2012.06.008

Lachaux, J.-P., Jerbi, K., Bertrand, O., Minotti, L., Hoffmann, D., Schoendorff, B., Kahane, P., 2007.
A Blueprint for Real-Time Functional Mapping via Human Intracranial Recordings. PLoS ONE
2, e1094. https://doi.org/10.1371/journal.pone.0001094

Lilly, J.M., Olhede, S.C., 2010. On the Analytic Wavelet Transform. IEEE Trans. Inform. Theory 56, 4135–4156. https://doi.org/10.1109/TIT.2010.2050935

Llorens, A., Trébuchon, A., Liégeois-Chauvel, C., Alario, F.-X., 2011. Intra-cranial recordings of brain activity during language production. Frontiers in psychology 2, 375.
https://doi.org/10.3389/fpsyg.2011.00375

Lopes da Silva, F., Van Rotterdam, A., 2005. Biophysical aspects of EEG and magnetoencephalogram generation, in: Electroencephalography: Basic Principles, Clinical Applications, and Related Fields. pp. 107–124.

Magnotti, J.F., 2020. RAVE: Comprehensive open-source software for reproducible analysis and visualization of intracranial EEG data 8.

Mallat, S., 2009. A Wavelet Tour of Signal Processing.

Maris, E., Oostenveld, R., 2007. Nonparametric statistical testing of EEG- and MEG-data , & 164, 177–190. https://doi.org/10.1016/j.jneumeth.2007.03.024

Medina Villalon, S., Paz, R., Roehri, N., Lagarde, S., Pizzo, F., Colombet, B., Bartolomei, F., Carron, R., Bénar, C.G., 2018. EpiTools, A software suite for presurgical brain mapping in epilepsy: Intracerebral EEG. Journal of Neuroscience Methods 303, 7–15.
https://doi.org/10.1016/j.jneumeth.2018.03.018

Mercier, M.R., Bickel, S., Megevand, P., Groppe, D.M., Schroeder, C.E., Mehta, A.D., Lado, F.A., 2017. Evaluation of cortical local field potential diffusion in stereotactic electro-encephalography recordings: A glimpse on white matter signal. NeuroImage 147, 219-232.
https://doi.org/10.1016/j.neuroimage.2016.08.037

Miller, K.J., denNijs, M., Shenoy, P., Miller, J.W., Rao, R.P.N., Ojemann, J.G., 2007. Real-time functional brain mapping using electrocorticography. NeuroImage 37, 504-507.
https://doi.org/10.1016/j.neuroimage.2007.05.029

Milstein, J., Mormann, F., Fried, I., Koch, C., 2009. Neuronal Shot Noise and Brownian 1/f2 Behavior in the Local Field Potential 4, 1-5. https://doi.org/10.1371/journal.pone.0004338

Nakai, Y., Jeong, J., Brown, E.C., Rothermel, R., Kojima, K., Kambara, T., Shah, A., Mittal, S., Sood, S., Asano, E., 2017. Three- and four-dimensional mapping of speech and language in patients with epilepsy. Brain 140, 1351-1370. https://doi.org/10.1093/brain/awx051

Oostenveld, R., Fries, P., Maris, E., Schoffelen, J.-M., 2011. FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Computational intelligence and neuroscience 2011, 156869. https://doi.org/10.1155/2011/156869

Piai, V., Dahlslätt, K., Maris, E., 2015a. Statistically comparing EEG/MEG waveforms through successive significant univariate tests: How bad can it be?: Comparing ERPs with successive univariate tests. Psychophysiol 52, 440-443. https://doi.org/10.1111/psyp.12335

Piai, V., Riès, S.K., Knight, R.T., 2015b. The electrophysiology of language production: what could be improved. Front. Psychol. 5. https://doi.org/10.3389/fpsyg.2014.01560

Riès, S.K., Dhillon, R.K., Clarke, A., King-Stephens, D., Laxer, K.D., Weber, P.B., Kuperman, R.A., Auguste, K.I., Brunner, P., Schalk, G., Lin, J.J., Parvizi, J., Crone, N.E., Dronkers, N.F., Knight, R.T., 2017. Spatiotemporal dynamics of word retrieval in speech production revealed by cortical high-frequency band activity. Proceedings of the National Academy of Sciences 114, E4530-E4538. https://doi.org/10.1073/pnas.1620669114

Roehri, N., Lina, J.-M., Mosher, J.C., Bartolomei, F., Benar, C.-G., 2016. Time-Frequency Strategies for Increasing High-Frequency Oscillation Detectability in Intracerebral EEG. IEEE Trans. Biomed. Eng. 63, 2595-2606. https://doi.org/10.1109/TBME.2016.2556425

Sahin, N.T., Pinker, S., Cash, S.S., Schomer, D., Halgren, E., 2009. Sequential processing of lexical, grammatical, and phonological information within Broca's area. Science (New York, N.Y.) 326, 445-9. https://doi.org/10.1126/science.1174481

Schalk, G., Leuthardt, E.C., Brunner, P., Ojemann, J.G., Gerhardt, L.A., Wolpaw, J.R., 2008. Real-time detection of event-related brain activity. NeuroImage 43, 245-249.
https://doi.org/10.1016/j.neuroimage.2008.07.037

Tadel, F., Baillet, S., Mosher, J.C., Pantazis, D., Leahy, R.M., 2011. Brainstorm: a user-friendly
application for MEG/EEG analysis. Computational intelligence and neuroscience 2011, 879716.
https://doi.org/10.1155/2011/879716

Tallon-Baudry, C., Bertrand, O., 1999. Oscillatory gamma activity in humans and its role in object representation. *Trends in cognitive sciences* 3, 151–162.

Trébuchon, A., Liegeois-Chauvel, C., Martinez, J.G., 2020. Contributions of Electrophysiology for Identifying Cortical Language Systems in Patients with Epilepsy. *Epilepsy & Behavior* 33.

Vidal, J.R., Freyermuth, S., Jerbi, K., Hamame, C.M., Ossandon, T., Bertrand, O., Minotti, L., Kahane, P., Berthoz, A., Lachaux, J.-P., 2012. Long-Distance Amplitude Correlations in the High Gamma Band Reveal Segregation and Integration within the Reading Network. *Journal of Neuroscience* 32, 6421–6434. https://doi.org/10.1523/JNEUROSCI.4363-11.2012

Vidal, J.R., Ossandón, T., Jerbi, K., Dalal, S.S., Minotti, L., Ryvlin, P., Kahane, P., Lachaux, J.-P., 2010. Category-Specific Visual Responses: An Intracranial Study Comparing Gamma, Beta, Alpha, and ERP Response Selectivity. *Frontiers in human neuroscience* 4, 195.
https://doi.org/10.3389/fnhum.2010.00195