# Signal processing and computational modeling for interpretation of SEEG-recorded interictal epileptiform discharges in epileptogenic and non-epileptogenic zones

Elif Koksal-Ersoz, Remo Lazazzera, Maxime Yochum, Isabelle Merlet, Julia

Makhalova, Borja Mercadal, Roser Sanchez-Todo, Giulio Ruffini, Fabrice

Bartolomei, Pascal Benquet, et al.

► To cite this version:

Elif Koksal-Ersoz, Remo Lazazzera, Maxime Yochum, Isabelle Merlet, Julia Makhalova, et al. Signal processing and computational modeling for interpretation of SEEG-recorded interictal epileptiform discharges in epileptogenic and non-epileptogenic zones. Journal of Neural Engineering, 2022, 19 (5), pp.055005. 10.1088/1741-2552/ac8fb4. hal-03796310

HAL Id: hal-03796310

https://hal.science/hal-03796310v1

Submitted on 4 Oct 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 CC BY 4.0 - Attribution - International License
PAPER • OPEN ACCESS

# Signal processing and computational modeling for interpretation of SEEG-recorded interictal epileptiform discharges in epileptogenic and non-epileptogenic zones

To cite this article: Elif Köksal-Ersöz et al 2022 *J. Neural Eng.* **19** 055005

View the [article online](#) for updates and enhancements.

You may also like

- [Collisional effect on the time evolution of ion energy distributions outside the sheath during the afterglow of pulsed inductively coupled plasmas](https://doi.org/10.1063/5.0059416)
J B Lee, H Y Chang and S H Seo

- [Ion energy and angular distributions in planar Ar/O₂ inductively coupled plasmas: hybrid simulation and experimental validation](https://doi.org/10.1063/5.0071371)
Yu-Ru Zhang, Zhuan-Zhuan Zhao, Chan Xue et al.

- [Measurements of ion energy distributions in a dual-frequency capacitively coupled plasma for Ar/O₂ discharges](https://doi.org/10.1063/5.0071372)
Jia Liu, Quan-Zhi Zhang, Yong-Xin Liu et al.

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

This content was downloaded from IP address 193.48.96.13 on 04/10/2022 at 12:11
<figure><img src="image_3.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

30 April 2022

REVISED

26 July 2022

ACCEPTED FOR PUBLICATION

6 September 2022

PUBLISHED

19 September 2022

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

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

PAPER

# Signal processing and computational modeling for interpretation of SEEG-recorded interictal epileptiform discharges in epileptogenic and non-epileptogenic zones

Elif Köksal-Ersöz¹,⁴,* ID, Remo Lazazzera¹,⁴ ID, Maxime Yochum¹ ID, Isabelle Merlet¹ ID, Julia Makhalova² ID, Borja Mercadal¹ ID, Roser Sanchez-Todo³ ID, Giulio Ruffini³ ID, Fabrice Bartolomei² ID, Pascal Benquet¹,⁴ ID and Fabrice Wendling¹,⁴,* ID

¹ Université de Rennes, INSERM, LTSI—UMR 1099, F-35000 Rennes, France

² Assistance Publique—Hôpitaux de Marseille, Service d'Epileptologie et de Rythmologie Cérébrale, Hôpital La Timone, Marseille, France

³ Neuroelectrics, Av. Tibidabo 47b, 08035 Barcelona, Spain

⁴ Equal contribution.

* Authors to whom any correspondence should be addressed.

E-mail: elif.koksal-ersoz@inserm.fr and fabrice.wendling@inserm.fr

**Keywords:** partial epilepsy, interictal epileptiform discharges, identification of epileptogenic zone, stereo-electroencephalography, neural mass models

## Abstract

**Objective.** In partial epilepsies, interictal epileptiform discharges (IEDs) are paroxysmal events observed in epileptogenic zone (EZ) and non-epileptogenic zone (NEZ). IEDs' generation and recurrence are subject to different hypotheses: they appear through glutamatergic and gamma-aminobutyric acidergic (GABAergic) processes; they may trigger seizures or prevent seizure propagation. This paper focuses on a specific class of IEDs, spike-waves (SWs), characterized by a short-duration spike followed by a longer duration wave, both of the same polarity. Signal analysis and neurophysiological mathematical models are used to interpret puzzling IED generation. **Approach.** Interictal activity was recorded by intracranial stereo-electroencephalography (SEEG) electrodes in five different patients. SEEG experts identified the epileptic and non-epileptic zones in which IEDs were detected. After quantifying spatial and temporal features of the detected IEDs, the most significant features for classifying epileptic and non-epileptic zones were determined. A neurophysiologically-plausible mathematical model was then introduced to simulate the IEDs and understand the underlying differences observed in epileptic and non-epileptic zone IEDs. **Main results.** Two classes of SWs were identified according to subtle differences in morphology and timing of the spike and wave component. Results showed that type-1 SWs were generated in epileptogenic regions also involved at seizure onset, while type-2 SWs were produced in the propagation or non-involved areas. The modeling study indicated that synaptic kinetics, cortical organization, and network interactions determined the morphology of the simulated SEEG signals. Modeling results suggested that the IED morphologies were linked to the degree of preserved inhibition. **Significance.** This work contributes to the understanding of different mechanisms generating IEDs in epileptic networks. The combination of signal analysis and computational models provides an efficient framework for exploring IEDs in partial epilepsies and classifying EZ and NEZ.

## 1. Introduction

During the pre-surgical evaluation of patients with drug-resistant epilepsy, stereo-electroencephalographic (SEEG, depth electrodes) recordings are performed

to identify regions involved in epileptic seizures. Identification of epileptogenic zones (EZs) and non-epileptogenic zones (NEZs, also referred to as propagation, irritative or non-involved zones) is usually based on pre-ictal and ictal patterns and ictal
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Example SEEG interictal activity recorded in a patient. The visualized SEEG recording is reduced to a small subset of bipolar channels for simplicity. In practice, such SEEG recordings are performed on a more significant number of recording channels (up to 128–256). The red and green markers correspond to spike-waves in the epileptogenic zone (EZ) and non-epileptogenic zone (NEZ). Electrodes are located as follows. H: Heschl gyrus; T: superior temporal gyrus (anterior part); A: amygdala (mesial contacts), middle temporal gyrus (lateral contacts); B: anterior hippocampus (mesial contacts), middle temporal gyrus (lateral contacts); C: posterior hippocampus (mesial contacts), middle temporal gyrus (lateral contacts); TP: temporal pole; TB: temporo-basal region; OF: orbito-frontal region.</figcaption></figure>

propagation. This categorization is used for deciding
the best therapeutic strategy.

Numerous studies have been dedicated to inter-
ictal paroxysmal events in humans [1]. Interictal
epileptiform discharges (IEDs) observed in SEEG
recordings demonstrate a wide range of morpho-
logy, ranging from simple monophasic epileptic
spikes to more complex multiphasic transient events.
Generation of IEDs in partial epilepsies is com-
monly ascribed to enhanced excitatory interactions
within glutamatergic neural networks. However, both
human data and animal models have supported the
view that gamma-aminobutyric acidergic (GABAer-
gic) signaling does play a significant role [2, 3]. The
presence of heterogeneous neural firing patterns dur-
ing IEDs in epileptic and non-epileptic zones sug-
gests a multitude of agents taking part in IED genera-
tion [4], including paroxysmal depolarization shifts
(PDSs) [5, 6] and depolarizing GABAergic actions
[7, 8]. So far, mechanisms leading to the genera-
tion and propagation of SEEG-recorded IEDs remain
unclear, and their clinical value for predicting the epi-
leptogenicity level of brain regions, is uncertain.

This paper focuses on a specific IED class called
biphasic epileptic spikes, typically observed in SEEG
recordings. These events, often referred to as spike-
waves (SWs), are characterized by a short-duration
spike followed by a longer duration wave, both with
the same polarity [9]. Figure 1 provides an example
of interictal activity recorded by SEEG electrodes
from a patient suffering from temporal lope epi-
lepsy (TLE). As part of a general clinical procedure,
the electrode implantation is specific to the patient

and it depends on the EZ location provided by non-
invasive methods: e.g. electroencephalogram (EEG),
MRI, and computed tomography scan. Figure 1 dis-
plays the selected bipolar channels where the IEDs
are less affected by volume conduction. We recall
that a bipolar channel captures the local potential
(defined as the electric potential difference between
adjacent contacts of an implanted electrode). After
this initial selection, SEEG experts labeled the bipolar
channels located in the EZ and NEZ. The red and
green markers correspond to the SWs in the EZ and
NEZ regions, respectively. By visual inspection, we
remark that spike and wave amplitudes are higher
in the EZ than in the NEZ with narrower spikes
and wider waves. These findings are supported by
literature by Alarçon et al [10] and more recently
by Serafini [11], both investigating the origin and
propagation of IEDs in the acute electrocorticogram.
In addition to these morphological differences, SWs
are more frequent in the EZ, and they can appear as
clusters and be preceded and/or followed by spiking
activity.

In the present manuscript, biphasic IEDs are
analyzed in-depth and the recorded IEDs are clas-
sified according to their morphologies, durations,
and generation sites. Their role in interictal epilep-
togenic networks is discussed accordingly to identify
the patient’s EZ versus NEZs. The discriminative role
of IEDs’ spatio-temporal features in the assessment
of the EZ and NEZ is compared to SEEG experts’
analysis. Subsequently, the physiological hypotheses
behind IEDs morphology are formulated using lam-
inar neural mass models (NMMs).
## 2. Methods

### 2.1. Patients and SEEG recordings

Five patients (labeled from P1 to P5) undergoing pre-surgical evaluation of drug-resistant partial epilepsy were selected for the study. Their clinical, anatomical, and electrophysiological data showed TLE in which parahippocampal and neocortical structures were involved. The group of patients in whom we analyzed SWs was homogeneous in terms of underlying pathology. The selection was based on precise criteria. First, patients were MRI-negative showing no focal cortical dysplasia. Second, very frequent SWs were recorded at electrode contacts located in mesial temporal lobe structures as well as in temporal neocortical structures.

SEEG recordings were performed using intracerebral multiple contact electrodes (10–15 contacts, length: 2 mm, diameter: 0.8 mm, 1.5 mm apart) placed following Talairach's stereotactic method [12, 13]. The positioning of electrodes was determined based on the available non-invasive information and clinical hypotheses about the localization of the EZ. The SEEG signals were recorded on a Deltamed system with a maximum of 128 channels. They were sampled at 256 Hz using no digital filter. For each patient, a mean of five SEEG recordings were analyzed (see table 1 for details). For each patient, interictal and ictal (where at least one seizure was recorded) periods were present. Seizures recorded during presurgical video-EEG monitoring were relatively reproducible although drug withdrawal had an impact on seizure patterns, as expected. The SEEG experts used the ictal periods to assess the epileptogenicity level of the selected bipolar channels. The interictal periods were exploited for the SWs detection and the SWs features extraction in order to affiliate the bipolar channel to the EZ versus the NEZ.

### 2.2. SEEG signals analysis

#### 2.2.1. SWs detection

A bipolar montage was adopted for the SEEG analysis, and prior to the investigation, a notch finite impulse response (FIR) filter was applied to remove the 50 Hz power supply component from signals. Among the IEDs detection approaches [14], the one proposed in [15] was implemented in the present work. This approach was adopted because it compares the input signal to a family of wavelets whose shape resembles the SWs morphology (complex Mexican hat wavelet). This comparison allows for enhancing the detection of IEDs in the records. In more detail, this method uses the mean value of the squared modulus of a wavelet filter banks output to increase the signal-to-noise ratio, as in [16]. The obtained signal was then used as input to the Page-Hinkley algorithm [17, 18]. This method builds a monotonically decreasing statistic, sample by sample from the input signal. In correspondence to an SWs, the behavior of the statistic

changes because the signal is above its mean: in fact, the metric starts to increase. It is at this point that the algorithm decides whether a spike is present, by comparing the increase of the statistic with a pre-defined threshold. The algorithm parameters were adjusted by visual inspection to minimize false-negative and false-positive detections.

#### 2.2.2. Bipolar channels selection

Generally speaking, SEEG recordings make use of multiple electrodes (usually between five and ten), and each electrode contains up to fifteen contacts. Among all the bipolar channel recordings, the objective of this step is to retain only those pairs of electrodes contact revealing a large number of IEDs. Some contacts may record the same IEDs generators: to avoid this problem the same strategy was adopted for all five patients considered in the study. It consisted in discarding bipolar channels containing IEDs generated by the same generator (usually adjacent plots referring to similar brain structures) or generated by volume conduction, and retaining only bipolar channels with a maximum number of IEDs. Table 1 provides the average number of SWs detected per interictal record and patient. On average, 365 000 SWs per patient were detected and exploited in the first screening. It is worth mentioning that the SWs distribution across the bipolar channels was not uniform. For this reason, the first screening of bipolar channels was performed by selecting only those showing the most significant number of IEDs with respect to the adjacent SEEG contacts.

#### 2.2.3. EZ vs. NEZ bipolar channels discrimination

A team of SEEG experts classified the selected bipolar channels of SEEG recordings performed in each patient as being located in the epileptogenic (EZ) [19] versus propagation or not involved zone (NEZ), by visual inspection. Ictal data was used to determine the degree of 'epileptogenicity' of a given region and its subsequent contribution to the EZ. The experts considered three criteria: the capability of a given region to generate high-frequency oscillations (rapid discharges) at seizure onset [18, 20–22], the delay of involvement of the region to the seizure onset and the seizure-related signal amplitude. SEEG experts worked independently and afterward, an agreement was reached to choose two bipolar channels per patient: one located in the EZ and another one in the NEZ.

Table 2 summarizes the selection of bipolar channels by SEEG experts. Except for A1–A2 and B1–B2, all channels are located in neocortical structures. Table 2 also provides the number of SWs detected per patient and bipolar channel. It is possible to observe the occurrence variability of interictal events for each patient and a higher number of SWs in the EZ than in the NEZ (except for P3).
<table><caption>**Table 1.** Information of ictal and interictal recordings per patient.</caption><thead><tr><th>Patient</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th></tr></thead><tbody><tr><td>ICTal records</td><td>2</td><td>4</td><td>1</td><td>4</td><td>1</td></tr><tr><td>Average ictal records duration (min)</td><td>25</td><td>38</td><td>60</td><td>22</td><td>30</td></tr><tr><td>Inter-ictal records</td><td>5</td><td>2</td><td>9</td><td>6</td><td>5</td></tr><tr><td>Average inter-ictal records duration (min)</td><td>52</td><td>60</td><td>56</td><td>53</td><td>60</td></tr><tr><td>SEEG contacts</td><td>93</td><td>116</td><td>95</td><td>116</td><td>62</td></tr><tr><td>Detected SWs per bipolar channel over all interictal records</td><td>3850</td><td>555</td><td>7103</td><td>3987</td><td>4269</td></tr><tr><td>Total detected SWs on interictal records</td><td>358 024</td><td>64 402</td><td>674 809</td><td>462 518</td><td>264 683</td></tr></tbody></table>

<table><caption>**Table 2.** SEEG experts selected one bipolar channel located in the EZ and one positioned in the NEZ for each patient. Electrode locations: TB, temporal basalis area; TP, temporal pole; A, amygdala complex; PH: parahippocampal gyrus; B, anterior part of Hippocampus. The number SWs detected per bipolar channel is reported together with the total of SWs detected per patient in both EZ and NEZ.</caption><thead><tr><th>Patient</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th></tr></thead><tbody><tr><td>Zone</td><td>EZ</td><td>n-EZ</td><td>EZ</td><td>n-EZ</td><td>EZ</td></tr><tr><td>Channel</td><td>TB1-TB2</td><td>TB9-TB10</td><td>TP2-TP3</td><td>TB5-TB6</td><td>A1-A2</td></tr><tr><td>SWs</td><td>4918</td><td>2386</td><td>2679</td><td>607</td><td>7821</td></tr><tr><td>SWs Sum</td><td></td><td>7304</td><td></td><td>3286</td><td>18 827</td></tr></tbody></table><table><thead><tr><th>Patient</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th></tr></thead><tbody><tr><td>Zone</td><td>EZ</td><td>n-EZ</td><td>EZ</td><td>n-EZ</td><td>EZ</td></tr><tr><td>Channel</td><td>TB1-TB2</td><td>TB9-TB10</td><td>TP2-TP3</td><td>TB5-TB6</td><td>A1-A2</td></tr><tr><td>SWs</td><td>4918</td><td>2386</td><td>2679</td><td>607</td><td>7821</td></tr><tr><td>SWs Sum</td><td></td><td>7304</td><td></td><td>3286</td><td>18 827</td></tr></tbody></table>

### 2.2.4. SWs mean waveform computation

The last step consisted of computing the mean wave-forms of the detected SWs per selected bipolar channel per record. This task aimed to identify qualitative morphological differences among SWs. For each record, segments of 1.5 s centered on the SWs detections were extracted and normalized using the z-score. A first SW mean waveform was obtained by averaging all channel detections. Then, each detection was aligned to the maximum of the cross-correlation function, computed pairwise between the detection and the preliminary average waveform. The final waveform was computed by averaging all the aligned detections. As shown in table 1 several thousands of SWs (3953 on average) were exploited to compute the mean SW waveform per bipolar channel.

By averaging waveforms of IEDs over time, we assumed that the morphology of these is stationary over time. We provided some evidence of this assumption by also visualizing a small subset of SWs (n = 30), randomly chosen across the time range for the same bipolar channels.

### 2.2.5. SWs features extraction and statistical tests

For each patient, we randomly selected 100 SWs at the EZ bipolar channel and 100 SWs at the NEZ bipolar channel to extract the SW features. If more than one interictal period was recorded in a patient, we randomly extracted the same amount of SWs from each period. Bipolar channels were visually inspected to avoid IEDs polarity inversion: when it was the case, the signal was inverted for SWs features extraction. An expert visually verified each SW in the set of 100 samples. False detections were discarded and replaced by the correct ones. For a selected SW, several features were extracted: spike and wave amplitude (calculated from the baseline depicted in red, in figure 2),

spike and wave full width at half maximum (FWHM, related to the sharpness of the spike and duration of the wave), spike peak to wave peak time delay (SW Delay), spike FWHM to wave FWHM time delay (FWHM Delay), the ratio between FWHM wave and FWHM spike, the ratio between Spike Amplitude and Wave Amplitude and the ratio between FWHM wave and FHWM Delay. The selected features are shown in figure 2 and enumerated in table 3.

Features were z-scored per patient: each feature value was normalized by subtracting its mean and dividing by its standard deviation, both computed for every patient. For example, a Spike-Wave Delay value equal to -0.9 indicates that the time delay between the spike and wave peaks is almost -1 standard deviation less with respect to the mean value of the delay variable computed on the patient data. This step was introduced to compare SWs features

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Extracted featured per SW: spike and wave amplitude, spike and wave full width at half maximum (FWHM), spike peak to wave peak time delay (SW Delay), spike FWHM to wave FWHM time delay (FWHM Delay), the ratio between FWHM wave and FWHM spike, the ratio between Spike Amplitude and Wave Amplitude and the ratio between FWHM wave and FHWM Delay.</figcaption></figure>
<table><caption>**Table 3.** Spike and wave extracted features.</caption><thead><tr><th>Features</th><th>Unit</th><th>Abbreviation</th></tr></thead><tbody><tr><td>Spike Amplitude</td><td>microVolt (μV)</td><td>Spike Amp</td></tr><tr><td>Wave Amplitude</td><td>microVolt (μV)</td><td>Wave Amp</td></tr><tr><td>Spike-Wave Delay</td><td>second (s)</td><td>SW Delay</td></tr><tr><td>Spike Full Width at Half Maximum</td><td>second (s)</td><td>FWHM Spike</td></tr><tr><td>Wave Full Width at Half Maximum</td><td>second (s)</td><td>FWHM Wave</td></tr><tr><td>Delay between FWHM Spike and FWHM Wave</td><td>second (s)</td><td>FWHM Delay</td></tr><tr><td>Ratio between Spike Amp and Wave Amp</td><td>—</td><td>Spike Amp/Wave Amp</td></tr><tr><td>Ratio between FWHM Wave and FWHM Spike</td><td>—</td><td>FWHM Wave/FWHM Spike</td></tr><tr><td>Ratio between FWHM Wave and FWHM Delay</td><td>—</td><td>FWHM Wave/FWHM Delay</td></tr></tbody></table>

extracted from different patients. Thanks to the z-score, we did not compare the absolute value of features but normalized deviations from their mean. This allowed considering patient-specific physiology. In other words, we were able, for example, to assess similarity in the deviation from the mean for the feature FWHM Spike, in the EZ and NEZ across different patients. The same comparison, but using features absolute values would have been more difficult because of inter-patient SWs morphology variability.

We used the Wilcoxon signed-rank test and the Kolmogorov-Smirnov test to compare the EZ versus the NEZ features distributions. The Kolmogorov-Smirnov test verifies the null hypothesis that the EZ (and NEZ) samples distribution follows a Gaussian law for a given feature. The Wilcoxon signed-rank test verifies the null hypothesis that two related paired samples come from the same distribution. A *p*-value less than 0.05 gives confidence of 95% to reject the null hypothesis. We tested if feature distributions followed a Gaussian law or not (Kolmogorov Smirnov test), and selected a nonparametric test (Wilcoxon signed-rank test) instead of Student's *t*-test in the latter case. Tests were performed in both ways: per patient and by merging all patients' SWs features data together.

### 2.2.6. SWs features clustering

We used k-means unsupervised clustering to discern if the SW morphology could be discriminative of the affiliation of SW to the EZ or NEZ. The clustering problem was performed on combined data from all patients and on data per patient to understand patient-specific features variations. The k-means clustering algorithm was applied with default parameters [23] by setting the number of clusters to two. In order to reduce the grouping dimensionality of the proposed features, we performed the clustering task by using tuples of three and two features.

Our interest was to establish if interictal activity can reveal concordant information on the EZ and NEZ to partition the data in a manner similar to the SEEG experts' analysis. In fact, the choice of EZ and NEZ was performed a priori on seizure data and then compared to the analysis performed on interictal data.

## 2.3. Computational model

NMMs are mathematical representations of the dynamics of the mean activity of synaptically connected neuronal subpopulations [24, 25]. The neural mass formalism assumes that these populations are almost homogeneous and 'synchronized', therefore, their behaviors can be reflected by average firing rates and membrane perturbations. NMMs have been extensively used to study both physiological [26] and pathological brain rhythms [27–29]. In the NMM considered in the present study, the activity of each subpopulation is given by a sigmoid 'pulse to wave' function that transforms the synaptic inputs *v* into a firing rate:

$$S(v) = \frac{5}{1 + e^{0.56(5 - v)}}. \quad (1)$$

A second-order differential equation relates the average presynaptic firing rate to post-synaptic potential (PSP):

$$\ddot{y} = \frac{W}{\tau_w} S(v) - \frac{2}{\tau_w} \dot{y} - \frac{1}{\tau_w^2} y, \quad (2)$$

where *W* represents the average synaptic gain, and $\tau_w$ is the synaptic time constant lumping the rise and decay times (assumed equal), and axonal delays. The time constant parameter $\tau_w$ reflects the kinetics of glutamatergic and GABAergic PSPs of cortical neurons. Given that an NMM formulates the 'synchronized' average temporal activity of a homogeneous population, single cell synaptic time constants can be conveyed (to some extent) to the model. On the other hand, parameters like average synaptic gain and connectivity coefficients are rather kept at an abstract level (they are 'lumped' parameters collecting various aspects not explicitly described by the model). Appendix D provides the system equations and parameters used in the manuscript.

The generation of realistic SEEG data to compare with real recordings requires combining the NMM framework with a physical formalism that can account for the location of the SEEG electrode contacts, the biophysical features of the volume conductor (tissue morphology and conductivity), and the physics of electrical current propagation deriving from Maxwell's equation. The integration of these
aspects is crucial for solving the SEEG-forward problem and hence simulating realistic SEEG signals. The present model includes a laminar architecture that represents the cortical layers of the human brain, following the framework in [30, 31] and also used in [32] for the analysis of seizure data. A laminar representation of a human cortical column of six layers with a physiological thickness [33] and uniform conductivity $\sigma = 0.3 \times 10^{-3} \text{ S mm}^{-1}$ [34] was considered. Given the anatomical and spatial characteristics of the pyramidal cells, we assumed that only the synaptic inputs on the pyramidal cells contribute to the local field potential (LFP) signal measured by a virtual SEEG electrode [35]. Synapses on the pyramidal cells were considered as point contacts as an approximation of the spatial distribution of the synaptic locations. They were placed along a one-dimensional fiber passing through cortical layers. Physiological constraints regarding the distribution of the synaptic contacts across the layers were loosened, and average synaptic locations were considered. Nevertheless, this simplification is appropriate given the level of the mesoscopic NMM formulation.

To a large extent, extra-cellular field potential (also known as LFP) is caused by the synaptic interactions between sub-populations of neurons and interneurons in a given neural population [35, 36]. Briefly, activation of an excitatory (inhibitory) synapse leads to a net inward (outward) flow of cations to the neuron, which depolarizes (hyperpolarizes) the membrane potential. Inflow (outflow) of positive charges into the neuron causes a negative (positive) current at the level of synaptic contact (so-called active sink (active source)). Meanwhile, due to the depolarizing (hyperpolarizing) membrane potential, further intra- and extra-cellular ion currents move along the membrane and a positive (negative) current will develop elsewhere along the neuronal element (so-called passive source (passive sink)). Then, there will be an extra-cellular current $I(t)$ flowing from the source to the sink. The field potential measured by a contact with reference $V_{\infty} = 0$ is then the difference between the source and the sink,

$$V(t) = \frac{I(t)}{4\pi\sigma} \left( \frac{1}{r_1} - \frac{1}{r_2} \right), \quad (3)$$

where $r_1$ and $r_2$ are the distance to the measurement point of the sink and of the source in a homogeneous infinite media, respectively. The amplitude $I(t)$ equals the post-synaptic current at the level of the synapses, and it is computed from the PSP given by the solution to equation (2),

$$I(t) = \eta y(t), \quad (4)$$

where $\eta = 10^{-3} \text{ S}$ is a conversion factor relating the PSP to post-synaptic current [31, 32]. In reality, $\eta$ is a lumped parameter that reflects neuronal population

characteristics, such as cell density and morphology, but these factors were ignored.

We did not model the electrode geometry but considered the voltage at the center of two contact points 2 mm away from each other. We assumed that the virtual SEEG electrode was parallel to the cortical column at a 10 mm distance from the considered neural mass. The system was simulated using the Euler-Murayama numerical method with a step size $dt = 10^{-4} \text{ s}$ and the same noise vector following a normal distribution.

## 3. Results

In the present section, results are reported following the methods presentation order.

### 3.1. SEEG analysis

#### 3.1.1. SWs mean waveform computation

The mean waveform is computed for each patient's EZ and NEZ bipolar channels. This representation is used to qualitatively understand morphological similarities and differences among the SWs and highlight the presence of discriminating characteristics to identify the EZ. In figures 3(a) and (b), the mean waveform of the patient P1's EZ and NEZ bipolar channels is represented in red; in blue are depicted all the SWs used for the mean waveform computation. The mean waveform was computed from an average of 6432 SWs per bipolar channel.

To show the IEDs morphology is stationary over time, figures 3(c) and (d) display in green the mean SWs mean waveform computed on a subset of IEDs ($n = 30$) randomly extracted over time for the same bipolar channels. The green signal is similar to the red one depicted in figures 3(a) and (b): this visually confirms that SWs morphology does not vary so much and that an average is appropriate.

The computation of the SWs mean waveform reveals qualitative differences between interictal events generated in the EZ with respect to those generated in the NEZ. In particular, the FWHM Spike is less pronounced in the EZ and the time delay between spike and wave is shorter. Also, FWHM Wave appears larger in the NEZ than in the EZ.

#### 3.1.2. SWs features distribution and statistical tests

For each patient, 200 SWs were used for features extraction: 100 to identify EZ SWs morphology and 100 to discover NEZ SWs characteristics. Features extracted from combined patients' records show different distributions between the EZ and NEZ. Figure 4 summarizes these findings in violin plots (the red area refers to the EZ while the green one refers to the NEZ). Dotted horizontal lines represent the quartiles of the distribution. These results indicate that spike amplitude is lower in the EZ than in the
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Detected spike-waves (SWs) and mean waveforms in patient P1. Detected SWs (blue) on the epileptogenic zone (a) and non-epileptogenic-zone (b) bipolar channels are used to compute the corresponding mean waveform (red). In (c) and (d) are reported the same plots in (a) and (b) but for a restrained number of SWs ($n = 30$) randomly selected over time: the green line is the mean SWs waveform computed on those 30 SWs and its shape is comparable to the red one computed respectively in (a) and (b). The plotted SWs were normalized, resulting in a dimensionless amplitude.</figcaption></figure>

NEZ, while the opposite is observed for wave amplitude. The SW Delay and FWHM Delay are shorter for the EZ SWs than the NEZ SWs. The FWHM Spike distribution demonstrates sharper spikes in the EZ SWs than in the NEZ SWs. Almost no difference exists in the FWHM wave feature in these two regions. The ratio features highlight a higher median value for FWHM Wave over FWHM Spike in the EZ, indicating that the EZ spikes have a sharper form than the NEZ spikes. The ratio of FWHM Wave over FWHM Delay reflects a smaller FWHM Delay in the EZ. Finally, the ratio Spike Amplitude/Wave Amplitude informs about a lower amplitude of spikes with respect to the Wave Amplitude in the EZ, but a higher Spike Amplitude with respect to the Wave Amplitude in the NEZ. This finding could be linked to the fact that a stronger inhibition effect could be found in the EZ due to the high excitability of the region. These

general features are considered for the SWs modeling. Inter-patient differences still exist, and those are reported in appendix A.

Table 4 reports the mean value of the features extracted from the SWs. In particular, the reported values correspond to the SWs located in the EZ and NEZ per patient; finally, the mean value of the features computed across patients in the EZ and NEZ groups is reported. From the table it is possible to draw some considerations: Wave Amplitude seems to be higher in the EZ (with respect to the NEZ) for all patients; SW Delay, FWHM Spike, FWHM Delay, and FWHM Wave resulted, instead, to be lower in the EZ for four patients over five.

The null hypothesis of the Kolmogorov Smirnov test (which assumes no difference between the observed and a normal distribution) is always rejected: p-value smaller than 0.05. The Wilcoxon
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Features distribution between epileptogenic zone (red) and non-epileptogenic zone (green). Dotted horizontal lines represent the quartiles of the distributions. For abbreviations refer to table 3. Data is normalized according to z-score; for details, refer to section 2.2.5.</figcaption></figure>

<table><caption>**Table 4.** Mean values of features computed per patient in the epileptogenic zone (EZ) and non-epileptogenic zone (NEZ). The lowest feature mean value, between the EZ and NEZ, is highlighted in blue, the highest in orange. For abbreviations see table 3.</caption><thead><tr><th rowspan="2">Feature</th><th colspan="2">P1</th><th colspan="2">P2</th><th colspan="2">P3</th><th colspan="2">P4</th><th colspan="2">P5</th><th colspan="2">All data</th></tr><tr><th>EZ</th><th>NEZ</th><th>EZ</th><th>NEZ</th><th>EZ</th><th>NEZ</th><th>EZ</th><th>NEZ</th><th>EZ</th><th>NEZ</th><th>EZ</th><th>NEZ</th></tr></thead><tbody><tr><td>Spike Amp</td><td>3.526</td><td>2.603</td><td>3.166</td><td>4.817</td><td>3.175</td><td>2.359</td><td>2.030</td><td>3.858</td><td>3.151</td><td>3.748</td><td>3.010</td><td>3.477</td></tr><tr><td>Wave Amp</td><td>2.857</td><td>2.132</td><td>2.386</td><td>2.114</td><td>2.911</td><td>2.895</td><td>3.410</td><td>2.509</td><td>3.209</td><td>2.850</td><td>2.955</td><td>2.500</td></tr><tr><td>SW Delay</td><td>0.157</td><td>0.216</td><td>0.284</td><td>0.278</td><td>0.261</td><td>0.261</td><td>0.094</td><td>0.174</td><td>0.089</td><td>0.182</td><td>0.177</td><td>0.222</td></tr><tr><td>FWHM Spike</td><td>0.015</td><td>0.043</td><td>0.045</td><td>0.061</td><td>0.032</td><td>0.035</td><td>0.025</td><td>0.023</td><td>0.020</td><td>0.035</td><td>0.027</td><td>0.040</td></tr><tr><td>FWHM Delay</td><td>0.060</td><td>0.096</td><td>0.144</td><td>0.134</td><td>0.118</td><td>0.120</td><td>0.045</td><td>0.089</td><td>0.033</td><td>0.080</td><td>0.080</td><td>0.104</td></tr><tr><td>FWHM Wave</td><td>0.199</td><td>0.250</td><td>0.311</td><td>0.238</td><td>0.253</td><td>0.265</td><td>0.114</td><td>0.144</td><td>0.146</td><td>0.201</td><td>0.204</td><td>0.219</td></tr><tr><td>FWHM W/FWHM S</td><td>15.161</td><td>6.602</td><td>7.136</td><td>3.912</td><td>8.556</td><td>9.470</td><td>5.292</td><td>6.460</td><td>8.195</td><td>5.791</td><td>8.868</td><td>6.447</td></tr><tr><td>Spike Amp/Wave Amp</td><td>1.440</td><td>1.422</td><td>1.636</td><td>2.452</td><td>1.160</td><td>0.879</td><td>0.609</td><td>1.867</td><td>1.106</td><td>1.485</td><td>1.190</td><td>1.621</td></tr><tr><td>FWHM W/FWHM D</td><td>3.598</td><td>2.781</td><td>2.349</td><td>2.004</td><td>2.485</td><td>2.486</td><td>2.641</td><td>1.977</td><td>5.103</td><td>2.532</td><td>3.235</td><td>2.356</td></tr></tbody></table>

signed-rank test, instead, shows that the EZ and the NEZ distribution differ for each feature, as the *p*-value is again smaller than 0.05.

### 3.1.3. *SWs features classification*

Overall, our results suggest that the EZ SWs present some subtle morphological differences with respect to those generated in a NEZ. We conducted a k-means classification using both two and three SWs extracted features from the selected patient recordings. The best results for the two features classification are accomplished by coupling the FWHM Delay and the FWHM Spike. The best three-element tuple candidate is composed using the following features:

FWHM Delay, the spike FWHM, and spike amplitude. Performances of these classifications are reported in figure 5 and table 5.

The complete list of binary classification performances obtained by combining different SWs features in couples or three-element tuples, are reported in appendix B.

The same k-means binary classification was conducted per patient using the same features couple and three-element tuple. Figure 6 and table 6 visually and quantitatively report the classification results performed per patient using two and three classification features. These results highlight inter-patient differences in the binary classification task.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Visual representation of a two (a) and three (b) features binary classification. Red elements refer to data classified by SEEG experts into the epileptogenic zone (EZ), and green ones belong instead to the non-epileptogenic zone (NEZ). Points represent correct classification data, while crosses refer to classification errors. For abbreviations refer to table 3. Data are normalized according to z-score; for details, refer to section 2.2.5.</figcaption></figure>

<table><caption>**Table 5.** Accuracy values were obtained by classifying two and three spike-wave features into two groups. For abbreviations see table 3.</caption><thead><tr><th>Feature 1</th><th>Feature 2</th><th>Feature 3</th><th>Classification accuracy (%)</th></tr></thead><tbody><tr><td>FWHM</td><td>FWHM</td><td></td><td>72.6</td></tr><tr><td>Delay</td><td>Spike</td><td></td><td></td></tr><tr><td>FWHM</td><td>FWHM</td><td>Spike Amp</td><td>73.2</td></tr><tr><td>Delay</td><td>Spike</td><td></td><td></td></tr></tbody></table>

The complete list of binary classification performances per patient, using different combinations of SWs features in couples or three-element tuples, is reported in appendix C.

## 3.2. Modeling IEDs

The SEEG signal analysis reveals four main differences between spike amplitude, FWHM Spike, SW Delay, and FWHM Delay. For the modeling framework, we focus on temporal features, such as FWHM Spike, SW Delay, and FWHM Delay. As given in table 2, identified EZs and NEZs in patients are in neocortical structures, except the EZs of P3 and P5. For this reason, we consider the six-layered cortical structure of the human cortex [33]. In the following, first, a minimal model of a cortical column for yielding IEDs with an SW morphology will be introduced. Then a network of two unidirectionally coupled columns, one for EZ and one for NEZ, will be considered to understand the mechanisms leading to the temporal differences between the IEDs in each region.

### 3.2.1. A laminar NMM for SW IEDs

We consider the same neural mass modeling approach as described in [27] for the hippocampus.

In this study, the model is adapted to the layered structure of the neocortex and endowed with the physical properties of conductive media [30–32]. The model comprises two populations of pyramidal cells (PYR and PYR'), and populations of somatostatin expressing (SST+) and parvalbumin expressing (PV+) GABAergic interneurons. In this model, the PYR' population accounts for the collateral excitation among PYR neurons. The PV+ and SST+ interneurons are locally excited by the PYR, and in turn, they mediate fast GABAergic and slow GABAergic inhibitory PSPs (IPSPs), respectively. Finally, the PV+ interneurons are inhibited by the SST+ interneurons.

Figure 7(a) demonstrates the laminar NMM with representative cell bodies and synaptic contacts in a cortical column. The PYR population expands along the column from layer 5 to layer 1. The synaptic inputs on the PYR population and their return currents contribute to the LFP recorded by two virtual SEEG electrode contacts, E1 and E2. We consider two synaptic locations on the PYR population: layer 1 and layer 5. The connections in layer 1 represent the synapses on the apical dendrites. The connections in layer 5 ensemble the synapses on the basal areas (dendritic and perisomatic). The apical synaptic inputs generate basal return currents, and basal synaptic inputs generate apical return currents. In this minimal model, the PYR' population and the PV+ interneurons target the basal area. The SST+ interneurons target both basal and apical dendrites of the pyramidal cells [37]. The model has only one population of SST+ interneurons; however, the PSP kinetics depends on the targeted cell region [38, 39], with the apical dendritic kinetics being slower than the basal dendritic kinetics. Figure 7(b) displays an
<figure><img src="image_10.png" /><figcaption>**Table 6.** Accuracy values are obtained by classifying two and three spike-wave features into two groups per patient. For abbreviations see table 3.</figcaption></figure>

<table><thead><tr><th rowspan="2">Feature 1</th><th rowspan="2">Feature 2</th><th rowspan="2">Feature 3</th><th colspan="6">Classification accuracy (%)</th></tr><tr><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th><th>Tot</th></tr></thead><tbody><tr><td>FWHM Delay</td><td>FWHM Spike</td><td></td><td>85.5</td><td>80</td><td>59.5</td><td>55.5</td><td>96.0</td><td>75.3</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td>Spike Amp</td><td>86.5</td><td>91.5</td><td>59.5</td><td>79.5</td><td>93.5</td><td>82.1</td></tr></tbody></table>

example time series with IEDs and figure 7(c) their mean waveform. PSPs on the PYR subpopulation in figure 7(d) show that the dynamics of the basal SST+ mediated IPSPs (IPSP$_{SST-B}$) are faster than the apical SST+ mediated IPSPs (IPSP$_{SST-A}$).

The generation of interictal spiking is a stochastic process, that is, the system is far from a periodic activity in the parameter space (which refers to an ictal period [27, 28]). Nevertheless, stochastic perturbations can trigger spiking behavior because the system is hyperexcitable due to an imbalanced level of excitation/inhibition. There are many different parameter combinations for obtaining an imbalanced level of excitation/inhibition in the model, yet we chose to follow a neurobiologically relevant strategy as much as the mesoscopic level of the model allows. For instance, we account for the reduction of the dendritic inhibition due to the loss of SST+ interneurons [40, 41], and depolarizing GABA due

to chloride transporter down regulation [42, 43]. Both are reflected by relatively small values for the synaptic amplitudes of the IPSPs mediated by the SST+ and PV+ interneurons. Moreover, an imbalanced excitation/inhibition level is also induced in the model parameters by increasing the glutamergic PSP amplitude, as shown in experimental animal models [44, 45]. Below, we concentrate on the impact of the synaptic dynamics (synaptic time constants and amplitudes) on the waveforms, which we see as the key parameters controlling the morphology of the IEDs. Yet, these parameters could lead to different waveforms under different parameter sets.

We explore mean waveform morphology against variations synaptic time constants and amplitudes by keeping the same spiking frequency. We consider the mean waveform in figure 7(c) as a reference for the comparison. Figures 8(a) and (b) show the impact of the basal and apical SST+ mediated IPSP amplitudes.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Laminar neural mass model for spike-wave generation and example waveforms. (a) Representative positions of the cell bodies and synaptic contacts in a cortical column of six layers. A virtual SEEG electrode with contacts E1 and E2 is positioned 10 mm away from the one-dimensional dendrite on the left. (b) An example time series of the local field potential at the level of contacts E1 and E2 in a bipolar montage. (c) Mean waveform of the interictal epileptiform discharge in panel (b). (d) Post synaptic potentials on the PYR subpopulation mediated by the excitatory inputs (EPSP), basal SST+ projection (IPSP<sub>SST-B</sub>), apical SST+ projection (IPSP<sub>SST-A</sub>) and basal PV+ projection (IPSP<sub>PV</sub>). Parameters are given in table D1.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Effect of the parameter variations on the mean waveforms. The reference mean waveform is shown in figure 7(c). (a), (b) Impact of the basal and apical inhibitory postsynaptic potential (IPSP) amplitudes. Stronger basal IPSP (a) gives a small spike with a prominent negative inflection. Stronger apical IPSP (b) gives a prominent positive wave with a small spike. (c), (d) Impact of the apical IPSP synaptic timescale. The slow apical IPSP (c) delays the wave and increases the wave half width. The fast apical IPSP (d) shortens the spike-wave delay and obscures the negative inflection. (e), (f) Impact of the basal IPSP synaptic timescale. The fast basal IPSP (e) narrows the negative inflection and advances the wave. The slow basal IPSP (f) widens the spike, delays, and downsizes the wave. (g), (h) The impact of the excitatory postsynaptic potential (EPSP) synaptic timescale. The fast EPSP (h) shortens the spike duration compared with slow EPSP (g). Features color codes in figure 2 are respected. Parameters are given in appendix D table D1.</figcaption></figure>

With stronger basal IPSP amplitude, the negative inflection dominates the morphology (figure 8(a)). In contrast, a stronger apical IPSP amplitude produces a dominant positive wave that can mask the initial spike (figure 8(b)). For the SW-type waveforms obtained under this variation, the delay between the spike and wave does not vary significantly. This observation indicates that the balance between the amplitudes of

the basal and apical SST+ mediated IPSPs determines the IEDs morphology.

Figures 8(c)–(f) show the impact of the basal and apical SST+ mediated IPSP kinetics on the mean waveform. Eventually, slower apical IPSPs delay the positive wave and increase the wave duration (figure 8(c)), while faster apical IPSPs advance and shorten the wave component (figure 8(d)). The
wave component masks the sharp negative inflection following the positive spike. As for the basal IPSP kinetics, the FHMW Spike, SW Delay, and FHMW Delay decrease for faster basal IPSPs (figure 8(e)). As the basal IPSP kinetics slows down and approaches the apical IPSP kinetics, the negative inflection and positive wave become less visible, and the spike widens. Finally, under the variation of the excitatory postsynaptic potential (EPSP) timescales, spikes narrow down for faster synapses (figure 8(g) vs. figure 8(h)), but the SW Delay does not change significantly.

The investigation above demonstrates the impact of the PSP time constants and the balance between the SST+ mediated basal and apical IPSPs on the SW morphology. For instance, strong apical IPSP masks the spikes, and increasing the synaptic time constant of the apical GABAergic currents increases the delay between peaks. Synaptic time constants can show some variability [46–50]. In addition, the LFP time profile can reflect the distribution/activation of synaptic channels or the level of synchronized synaptic inputs. While these factors can be responsible for the temporal differences of the EZ’s and NEZ’s IEDs, there might be other physiological factors that can account for the protective roles of the NEZ’s IEDs [2]; for instance, the prolonged activity of inhibitory neurons by slow glutamatergic inputs. We address this question in the next section.

### 3.2.2. A minimal network of EZ and NEZ

The main difference between an EZ and a NEZ is that the EZ is a hyperexcited region, where seizure is initiated, and frequent interictal spikes can be observed. Hyperexcitability in the EZ can be a consequence of many different factors, such as depolarizing GABA, loss of inhibitory cells, or malfunctioning inhibitory signaling. On the other hand, the excitation/inhibition balance is preserved in the NEZ. As a consequence, in the NEZ seizures are not initiated and interictal spiking is less frequent. Furthermore, as demonstrated in section 3.1, compared to the EZ SWs, the NEZ SWs have different spatio-temporal features: spikes are wider, and waves are well-delayed.

Here we consider a network of two unidirectionally coupled systems, one standing for an EZ and one for a NEZ (figure 9(a)). The system representing the EZ is in a hyperexcitable mode, i.e. static depolarizing GABA<fnref n="5" /> and enhanced glutamatergic

activity contributes to the generation of spontaneous IEDs. The excitation/inhibition balance is preserved in the system representing the NEZ by the inhibitory interneurons, which appease the irritative impact of the excitatory projections from the EZ system. Therefore, the NEZ system does not undergo spontaneous discharges but responds to the perturbations from the EZ system. A slow glutamatergic signaling, which could be triggered by activation of slow glutamatergic receptors, such as ionotropic NMDAR or post-synaptic mGluR, and an excitatory input are introduced into the NEZ system. In terms of synaptic contacts, the structure in figure 7(a) is adopted for the EZ and NEZ systems, plus a basal slow-glutamatergic and an apical excitatory input for the latter. We consider the same values for the synaptic time constants in two systems, with the slow glutamatergic EPSP being five-times slower than the (fast-)EPSP action. Activities of the EZ and NEZ systems are projected on two distinct virtual electrodes (two contacts on each), which are assumed to be measuring only one system.

Figures 9(b)–(e) show the time series of the EZ and NEZ interictal activity and the mean waveforms. The low frequency spiking activity in the NEZ (figure 9(d)) is due to an intrinsically less excited NEZ system and a moderate level of interaction between the two systems. In addition, long-lasting inhibitory activity reduces consecutive spiking (explained below). The mean waveform of the EZ IEDs has an SW morphology with an initial sharp spike followed by a slow wave. The mean waveform of the NEZ IEDs starts with a small negative early activity that reflects the apical excitatory input from the EZ system. It is followed by an SW-type morphology with wider spike and wave components with longer delays than the EZ-SW mean waveform. The mechanism behind these differences is the slow glutamatergic signaling in the NEZ system. Wider spikes reflect the prolonged glutamatergic activity meditated by the slow glutamatergic neurotransmitters. Similarly, the slow glutamatergic neurotransmitters targeting the SST+ interneurons sustain their activation. The sustained basal GABAergic activity increases the SW delay, and the sustained apical GABAergic activity prolongs the wave. The long-lasting inhibition also avoids the NEZ system undergoing spiking in response to incoming IEDs from the EZ system. In this sense, the long-lasting inhibition has a protective role in the NEZ. The difference in the PSP dynamics between the EZ and NEZ are visible in figures 9(f) and (g), where we can notice the prolonged inhibition by the slow glutamatergic drive.

<footnote n="5">⁵ Depolarizing GABA is generated by chloride accumulation in the neuron that changes the reversal potential for GABAA receptors (shifting from −80 mV to −50 mV). This chloride accumulation can occur in two different conditions. First, the neuron has a reduced number of KCC2 chloride transporter due to downregulation, but still enough to be able to extrude chloride outside the cell when the firing rate of GABAergic interneuron and consequent GABA release is moderated. However, when the firing rate of interneurons is high the chloride accumulates into the cell and progressively shifts the reversal potential from hyperpolarizing to</footnote>

<footnote>depolarizing. This mechanism is called activity-dependent depolarizing GABA [51]. Second, when neurons have a very low number of KCC2 because the majority of transporter have been downregulated, chloride cannot be extruded and GABA is always depolarizing [42]. This mechanism is called sustained or static depolarizing GABA.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 9.** A minimal network of the epileptogenic zone (EZ) and non-epileptogenic zone (NEZ) and example waveforms. (a) The EZ system (top) is unidirectionally coupled to the NEZ system (bottom) via an excitatory input. Dark red arrows in the NEZ model represent the slow-glutamatergic connections. A virtual electrode with contact E1 and E2 measures the activity of the EZ system. A virtual electrode with contact F1 and F2 measures the activity of the NEZ system. (b) Time series of the EZ system at the E1 and E2 contacts level in a bipolar montage. (c) The mean waveform of the interictal epileptogenic discharges in panel (b). (d) Time series of the NEZ system at the level of F1 and F2 contacts in a bipolar montage. (e) The mean waveform of the interictal epileptogenic discharges in panel (d). Features color codes in figure 2 are respected. (f) Post synaptic potentials on the PYR subpopulation in EZ mediated by excitatory inputs (EPSP), basal SST+ projection (IPSP<sub>SST-B</sub>), apical SST+ projection (IPSP<sub>SST-A</sub>) and basal PV+ projection (IPSP<sub>PV</sub>). (g) Post synaptic potentials on the PYR subpopulation in NEZ mediated by fast excitatory inputs (EPSP<sub>fast</sub>), basal SST+ projection (IPSP<sub>SST-B</sub>), apical SST+ projection (IPSP<sub>SST-A</sub>), basal PV+ projection (IPSP<sub>PV</sub>) and slow excitatory inputs (EPSP<sub>slow</sub>). Parameters are given in appendix D table D2.</figcaption></figure>

## 4. Discussion

The generation of IEDs in partial epilepsies is commonly ascribed to enhanced excitatory interactions within glutamatergic neural networks. However, both human data and animal models support the view that GABAergic signaling plays a significant role [2, 3]. This paper has studied a sample set of IEDs recorded by SEEG electrodes in a group of partial epilepsy patients. The recorded IEDs have been classified with respect to their morphologies, duration, and generation sites. Extracted features were then used to classify the epileptogenicity of the recording region. Then, a laminar NMM has been used to mimic the spatiotemporal features of the IEDs, which have been classified with respect to their generation site. Laminar contributions of the glutamatergic and GABAergic PSPs and the positions of the synaptic contacts have been considered to the simulated LFP signals. Slow glutamatergic signaling and long-lasting

inhibition have been marked as distinctive mechanisms in generating SW discharges in the EZ and NEZ models.

### 4.1. SWs features selection, distribution and clustering

Five patients with TLE were selected for SEEG recordings. The IEDs were detected on each bipolar channel from the interictal data, and SEEG experts identified one channel in the EZ and one in the NEZ for each patient. In total, 64 317 SWs have been used to compute the SW mean waveform per patient and per selected bipolar channel. Its qualitative analysis has revealed a sharper nature of spikes in the EZ and a longer SW delay in the NEZ. Among different possibilities, one can hypothesize that it can be due to the hypersynchronous activity of the pyramidal cells in EZ, which discharge all together simultaneously. In contrast, the pyramidal cells in NEZ can discharge asynchronously (in line with a balanced excitation/
inhibition ratio), which would give wider spikes and
delayed waves.

To identify the morphological differences in the
EZ and NEZ SWs, the features (table 3) were extrac-
ted from 100 SWs per bipolar channel. The randomly
constructed sets of 100 SWs were built with the intent
to maximize their internal heterogeneity by using all
interictal records available per patient and randomly
choosing the same quantity of SWs per recording.
Each set of SWs was visually verified by an expert to
avoid false detection. This process assured the correct-
ness of extracted features but also limited the num-
ber of samples in the set. Among the extracted fea-
tures, we have chosen two types of time delay: spike
peak to wave peak (SW Delay) and spike FWHM to
wave FWHM delay (FWHM Delay). This choice relies
on the fact that wave peaks were more susceptible to
noise (interfering with measurements) than the rapid
apparition of the wave next to the spike.

From the qualitative analysis of SWs mean wave-
forms and the quantitative results of unsupervised
classification (tables in appendices B and C), the fol-
lowing features have been revealed to be the most
physiologically significant to discriminate the EZ and
NEZ: Spike Amplitude, FWHM Delay and FWHM
Spike. Amid these three, the Spike Amplitude has
been discarded because this feature is also related
to the position of the electrodes with respect to the
sources. Intuitively, Spike Amplitude might be cor-
related with EZ. Instead, our results suggest that it
is not the case; for example, NEZ produces wider
spikes for P2, P4, and P5. Using only the FWHM
Delay and FWHM Spike as binary classification fea-
tures, we have reached an accuracy of 72.6% in separ-
ating all patients' SWs in EZ and NEZ; and 75.3% for
a patient-specific approach. It is worth noticing that
only using these two features, the classification accur-
acy was higher than 80% for three patients over five.
We have considered these features for the modeling
study. Features, like the ratio between FWHM wave
and FWHM spike, the ratio between Spike Amplitude
and Wave Amplitude and the ratio between FWHM
wave, and FHWM Delay, are a non-linear combina-
tion of primary features, and their use was exploited
to discover if they could lead to better classification
performances. Finally, these features have not repres-
ented a real advantage for classification.

Some comments about patients’ specific features
distribution (figure in appendix A) are worth men-
tioning. Particularly for patients P2 and P3, it is visu-
ally difficult to affirm that the EZ and NEZ distribu-
tions of some features are different. We can consider
two causes: the smaller number of data points (100
points per distribution) and the patient-specific mor-
phology. In particular, patient P2 feature distribu-
tions have shown that the FWHM Wave was bigger
in the EZ SWs than in the NEZ SWs. This finding can
be since the EZ presents frequent trains of successive

spikes (polyspikes), which can induce a longer wave
duration, usually linked to inhibition mechanisms
or ionic accumulation [2]. It was evident that the
NEZ bipolar electrodes presented a less spiky activ-
ity. Classification has shown poor discrimination per-
formances in patient P3. We have observed that the
SWs with opposite polarities were recorded in the EZ
through visual inspection of signals. The reason can
be that the EZ electrode was situated between two
epileptic sources. In this work, we have chosen only
one polarity type for the analysis, which decreased the
number of analyzed IEDs and could have led to poor
classification performances in patient P3. For com-
parison with literature, Serafini and Loeb [52] have
shown that on human SEEG, IEDs at the center of epi-
leptic foci exhibit a prevalent sharp wave while those
at the periphery exhibit a relative enhancement of the
slow-wave, which possibly corresponds to surround
inhibition. In a more recent study [11], Serafini has
defined ‘red-spikes zone’ as those generating seizures
and ‘green-spike zone’ as those generating only IEDs.
The author found that the slow waves are weakened
in the red-spikes zone, and their kinetics shows more
rapid decay in that area. From the SW features distri-
bution plots in appendix A, it is possible to appreciate
in four patients out of five a larger wave in the NEZ
and a sharper spike in the EZ.

Overall, our results have shown that the SWs mor-
phology can be predictive of the EZ versus the NEZ.
Nevertheless, the present work’s limitations consist of
the fact that among the large variability of IEDs, only
SWs have been considered. It is possible that discrim-
inating features could be present in other different
morphologies worth to be analyzed. The presented
results are consistent for patients with a specific dis-
order (TLE). A larger cohort of patients would con-
solidate the presented outcome. Among the detected
SWs, this paper has not explored the co-occurrence of
IEDs between the EZ and NEZ bipolar channels. We
believe this information is essential to understanding
brain dynamics and would allow assessing the direc-
tionality of information flow among brain areas.

4.2. Spatially distributed synaptic interaction for
the SW morphology

In order to model the IEDs measured by SEEG
we have used a neurophysiologically-plausible lam-
inar NMM. The synaptic interaction at different
layers with synaptic time constants accounting for
local dendritic delays has permitted to simulate real-
istic SW morphologies. While the simple model was
responding to the SW morphology, slow glutamater-
gic signaling was included to account for the tem-
poral differences between the IEDs in the EZs and
NEZs. Slow glutamatergic signaling has increased the
spike duration by targeting the PYR population and
sustained a long-lasting inhibition that increases the
SW delay and the wave duration of the NEZ IEDs.
In our model, the SST+ interneurons play a crucial role in generating IEDs through different mechanisms. First, as reported by previous computational studies [27, 28, 32, 53], they control the epi-leptogenicity of the system. Second, the spatial distribution of their axons through cortical layers generates IPSPs with different dipoles that shape the SW morphology in SEEG detected signals. Third, they mediate GABAergic signaling at different time scales due to the distance between the apical and basal synaptic contacts in the PYR population.

Unlike the previous models [27, 31, 32], the SST+ interneurons in our model target both superficial and deep layers of the cortical column with different kinetic properties. This structure is physiologically grounded. Indeed, the SST+ interneurons (Martinotti and non-Martinotti cells) group dendritic targeting cells that differ in morphology, electrophysiology, and connectivity [54]. Martinotti cells have both local axonal arbor where their somas are located and long ascending axons in layer 1. Non-Martinotti cells have local axons in layer 4 and layers 5/6 with unknown connectivity for the latter. Non-Martinotti cells are faster in generating action potentials than Martinotti cells. Interestingly, the SST+ interneurons have facilitating excitatory synapses regardless of their subtypes, and they can produce long-lasting spiking. In this stark contrast to other types of interneurons, a single high-frequency burst from an excitatory presynaptic cell can activate the SST+ interneurons and enable long-lasting feedback inhibition. Finally, recent studies have argued that the SST+ interneurons gate not only GABA<sub>A</sub>-mediated inhibition but also slow GABA<sub>B</sub>-mediated inhibition and can silence excitatory transmission [55].

SST+ interneurons provide lateral feedback inhibition, which can affect a wide area. Lateral inhibition mediated by Martinotti cells is facilitated by cholinergic inputs [56]. In the rodent hippocampus, balanced cholinergic and mGluR activity enhances and prolongs the activation of the SST+ interneurons [57]. Moreover, depolarization of Martinotti cells by group I mGluRs has been shown in the human neocortex [58], and this mechanism has been suggested as an antiepileptic pathway. Facilitating excitatory inputs and prolonged activation of the SST+ interneurons can have a crucial role in preserving the excitatory/inhibitory balance. We have not considered short-term plasticity, but it is plausible that facilitating excitatory synapses on the SST+ interneurons contribute to long-lasting GABAergic inhibition, hence, to SW generation.

PDS has been correlated with extra-cellular SW activity in animal models of epilepsy and in human laminar recordings [2, 5, 6]. Initialization of PDS has been correlated with excitatory synaptic currents that activate GABAergic responses. Ionic currents are also

essential contributors to the sequential occurrence of PDS events. Our model does not account for non-synaptic events but EPSP and IPSP. The chain of events in PSP dynamics leading to SWs, is in line with the ones of PDS generation. As for PDS synchronization within a population, it has been proposed to depend on both excitatory [2] and GABA-mediated depolarization [7, 59–61]. Depolarizing GABA in our simulations contributes to the hyperexcitability of the EZ and the spike amplitude of SWs in this region but not directly to their initialization. It would be interesting to investigate the GABA-mediated static depolarization in an extended model with recurrent interactions between the PV+ interneurons.

In terms of limitations, this study has assumed point contacts for synaptic interactions between different cell types. However, synaptic contacts are distributed all over the apical and basal dendrites of the principal cells. It is envisaged to use realistic compartmental models [62] to better model IEDs and to test different hypotheses. Besides, the contribution from the extra-cellular ionic concentrations or the activity of the glial cells to IEDs is not considered. These points will be considered in the future.

Another limitation of this study is modeling the cortical column with only one PYR population contributing to the LFP signal. While layer 5 pyramidal neurons generate the largest dipole in a cortical column, layer 2/3 or layer 6 pyramidal neurons might contribute to SEEG signals. Additionally, intra-laminar recordings have shown a laminar hierarchy in interictal spiking [5] and seizure propagation [63]. Other experimental and modeling studies have suggested laminar organization in SW generation [64, 65]. Eventually, models with different levels of cortical and subcortical organizations can generate richer IED patterns. For instance, a laminar based approach for simulating SW-type discharges in EEG recordings has been proposed in [66] where EEG output has been modeled as a weighted linear combination of the dipoles derived from the PSPs on pyramidal neuronal population. In this study, the author has shown that different EEG waveforms can be obtained by varying the dipole weights. Another important point in modeling SEEG is the effect of the position of the electrode contacts on the observed signal morphology, especially depending on the location of the electrode contacts being in the subcortical, white matter, or gray matter, which are ignored in the present manuscript. Nevertheless, we believe that the mechanism proposed in this work, that is apical and basal interactions generating opposite current dipoles leading to SW-type discharges, is general enough to be applied to structures other than the neocortex. Future works envisage investigating complex cortical interactions in different regions during interictal and ictal events including the role of electrode positions.
## 5. Conclusion

The relation between underlying pathology and the type of IEDs recorded in patients with partial epilepsy may be complex. This study has focused on a subset of IEDs, which are SWs discharges. IEDs have been detected from interictal data recorded in the EZ and NEZ, which were identified by the SEEG experts from the ictal data and fast onset activity. Some of the SW morphology features have been found to be discriminating of the EZ vs. the NEZ. In particular, SWs with a sharp spike and a small delay between spike-to-wave peaks are most likely found in bipolar channels that recorded fast onset activity. Our computational model has made a link between the morphology and the extracted information from the signal analysis and underlined the impact of inhibition. Our results suggest that interictal epileptogenic spikes can be indicative of the EZ and NEZ.

### Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.

### Acknowledgments

Conceptualization: E K-E, R L, P B and F W. Data curation: I M, J S, F B, P B and F W. Methodology:

E K-E, R L, B M, R S-T and G R. Funding accusation: G R, F B and F W. Formal analysis, research, validation, visualization, writing—original draft: E K-E and R L. Software: E K-E, R L, and M Y. Writing—review and editing: E K-E, R L, G R, P B and F W.

B M and R S-T work for Neuroelectrics Barcelona. G R is co-founder of Neuroelectrics.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 855109).

### Ethical statement

The SEEG recordings were carried out as part of normal clinical care of patients. All patients were informed that their data could be used for research purpose and they gave their written consent in the usual way.

## Appendix A. SW features distribution per patient

Feature distributions of EZ SWs (red) and NEZ SWs (green) are presented per patient. Dotted horizontal lines represent the quartiles of the distribution. For abbreviations refer to table 3. Data are normalized according to z-score; for details, refer to section 2.2.5.
(P1)

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

Spike Amp

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

Wave Amp

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

SW Time Delay

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

FWHM Delay

FWHM Spike

FWHM Wave

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

FWHM Wave
FWHM Spike

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

Spike Amp
Wave Amp

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

FWHM W
FWHM Delay

(P2)

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

Spike Amp

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

Wave Amp

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

SW Time Delay

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

FWHM Delay

FWHM Spike

FWHM Wave

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

FWHM Wave

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

FWHM Wave
FWHM Spike

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

Spike Amp
Wave Amp

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

FWHM W
FWHM Delay

(P3)

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

Spike Amp

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

Wave Amp

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

SW Time Delay

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

FWHM Delay

FWHM Spike

FWHM Wave

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

FWHM Wave
FWHM Spike

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

Spike Amp
Wave Amp

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

FWHM W
FWHM Delay

(P4)

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

Spike Amp

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

Wave Amp

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

SW Time Delay

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

FWHM Delay

FWHM Spike

FWHM Wave

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

FWHM Wave

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

FWHM Wave
FWHM Spike

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

Spike Amp
Wave Amp

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

FWHM W
FWHM Delay

((P5))

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

Spike Amp

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

Wave Amp

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

SW Time Delay

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

EZ

NEZ

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

FWHM Delay

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

FWHM Spike

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

FWHM Wave

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

FWHM Wave
FWHM Spike

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

Spike Amp
Wave Amp

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

FWHM W
FWHM Delay
## Appendix B. Classification results on all patients data

K-means binary classification accuracy using different SWs features combinations. For abbreviations refer to table 3.

<table><thead><tr><th>Feature 1</th><th>Feature 2</th><th>Feature 3</th><th>Classification accuracy (%)</th></tr></thead><tbody><tr><td>Wave Amp</td><td>FWHM Wave</td><td></td><td>54.4</td></tr><tr><td>FWHM Spike</td><td>FWHM Wave</td><td></td><td>58.1</td></tr><tr><td>Spike Amp</td><td>FWHM Wave</td><td></td><td>58.4</td></tr><tr><td>SW Delay</td><td>FWHM Wave</td><td></td><td>58.8</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td></td><td>60.9</td></tr><tr><td>FWHM Delay</td><td>FWHM Wave</td><td></td><td>61.2</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td></td><td>61.3</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td></td><td>63.6</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td></td><td>65.0</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td></td><td>65.7</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td></td><td>66.4</td></tr><tr><td>Spike Amp</td><td>FWHM Spike</td><td></td><td>67.8</td></tr><tr><td>Wave Amp</td><td>FWHM Spike</td><td></td><td>69.4</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td></td><td>72.6</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td></td><td>74.3</td></tr><tr><td>SW Delay</td><td>FWHM W/FWHM S</td><td></td><td>51.5</td></tr><tr><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td></td><td>56.6</td></tr><tr><td>Wave Amp</td><td>FWHM W/FWHM Delay</td><td></td><td>61.1</td></tr><tr><td>Spike Amp</td><td>FWHM W/FWHM Delay</td><td></td><td>61.4</td></tr><tr><td>Spike Amp/Wave Amp</td><td>FWHM W/FWHM Delay</td><td></td><td>61.7</td></tr><tr><td>Spike Amp</td><td>FWHM W/FWHM S</td><td></td><td>62.0</td></tr><tr><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td></td><td>62.7</td></tr><tr><td>Wave Amp</td><td>FWHM W/FWHM S</td><td></td><td>62.8</td></tr><tr><td>FWHM Delay</td><td>Spike Amp/Wave Amp</td><td></td><td>64.6</td></tr><tr><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td></td><td>66.4</td></tr><tr><td>SW Delay</td><td>Spike Amp/Wave Amp</td><td></td><td>66.8</td></tr><tr><td>SW Delay</td><td>FWHM W/FWHM Delay</td><td></td><td>67.7</td></tr><tr><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td></td><td>69.6</td></tr><tr><td>FWHM Spike</td><td>FWHM W/FWHM Delay</td><td></td><td>74.7</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM Wave</td><td>58.2</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM Wave</td><td>58.5</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM Wave</td><td>58.7</td></tr><tr><td>Wave Amp</td><td>FWHM Spike</td><td>FWHM Wave</td><td>61.1</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM Wave</td><td>61.8</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td>FWHM Wave</td><td>62.0</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM Wave</td><td>63.8</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>FWHM Wave</td><td>64.2</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM Delay</td><td>66.2</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM Delay</td><td>66.8</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM Delay</td><td>67.1</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td>FWHM Wave</td><td>67.1</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>SW Delay</td><td>67.5</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM Spike</td><td>67.9</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>FWHM Spike</td><td>68.1</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td>FWHM Wave</td><td>71.1</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td>FWHM Spike</td><td>72.6</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM Spike</td><td>73.2</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM Spike</td><td>73.5</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM Spike</td><td>75.0</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM W/FWHM S</td><td>50.4</td></tr><tr><td>SW Delay</td><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td>50.8</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM W/FWHM S</td><td>56.7</td></tr><tr><td>SW Delay</td><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td>58.4</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM W/FWHM Delay</td><td>63.2</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>64.2</td></tr></tbody></table>

(Continued.)
(Continued.)

<table><thead><tr><th>Feature 1</th><th>Feature 2</th><th>Feature 3</th><th>Classification accuracy (%)</th></tr></thead><tbody><tr><td>FWHM Delay</td><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td>64.3</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM W/FWHM S</td><td>64.8</td></tr><tr><td>FWHM Spike</td><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td>65.5</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM W/FWHM Delay</td><td>66.1</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>Spike Amp/Wave Amp</td><td>66.3</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>66.8</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>67.3</td></tr><tr><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td>67.9</td></tr><tr><td>Wave Amp</td><td>FWHM Spike</td><td>FWHM W/FWHM Delay</td><td>68.2</td></tr><tr><td>SW Delay</td><td>Spike Amp/Wave Amp</td><td>FWHM W/FWHM Delay</td><td>68.3</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM W/FWHM Delay</td><td>68.6</td></tr><tr><td>Spike Amp</td><td>FWHM Spike</td><td>FWHM W/FWHM Delay</td><td>69.3</td></tr><tr><td>FWHM Spike</td><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td>69.4</td></tr><tr><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td>FWHM W/FWHM Delay</td><td>70.7</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td>74.2</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td>FWHM W/FWHM Delay</td><td>74.7</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td>75.0</td></tr></tbody></table>

## Appendix C. Classification results per patient

K-means binary classification accuracy, per patient, using different SWs features combinations. For abbreviations refer to table 3.

<table><thead><tr><th>Feature 1</th><th>Feature 2</th><th>Feature 3</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th><th>Acc</th></tr></thead><tbody><tr><td>Wave Amp</td><td>FWHM Wave</td><td></td><td>53.0</td><td>68.0</td><td>53.0</td><td>66.0</td><td>62.5</td><td>60.5</td></tr><tr><td>SW Delay</td><td>FWHM Wave</td><td></td><td>66.0</td><td>66.0</td><td>51.0</td><td>65.5</td><td>70.5</td><td>63.8</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td></td><td>74.0</td><td>53.0</td><td>53.5</td><td>80.5</td><td>72.5</td><td>66.7</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td></td><td>71.5</td><td>55.5</td><td>52.0</td><td>71.0</td><td>90.5</td><td>68.1</td></tr><tr><td>FWHM Delay</td><td>FWHM Wave</td><td></td><td>74.5</td><td>72.0</td><td>53.0</td><td>64.0</td><td>81.5</td><td>69.0</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td></td><td>74.0</td><td>51.5</td><td>50.5</td><td>77.5</td><td>92.5</td><td>69.2</td></tr><tr><td>Spike Amp</td><td>FWHM Wave</td><td></td><td>65.5</td><td>86.0</td><td>60.5</td><td>77.0</td><td>59.0</td><td>69.6</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td></td><td>69.0</td><td>74.0</td><td>63.0</td><td>82.5</td><td>61.5</td><td>70.0</td></tr><tr><td>FWHM Spike</td><td>FWHM Wave</td><td></td><td>75.0</td><td>82.0</td><td>51.5</td><td>56.0</td><td>86.0</td><td>70.1</td></tr><tr><td>Wave Amp</td><td>FWHM Spike</td><td></td><td>81.5</td><td>69.5</td><td>51.5</td><td>69.0</td><td>86.5</td><td>71.6</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td></td><td>84.0</td><td>73.5</td><td>52.5</td><td>68.5</td><td>93.0</td><td>74.3</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td></td><td>74.0</td><td>78.0</td><td>63.0</td><td>77.5</td><td>80.5</td><td>74.6</td></tr><tr><td>Spike Amp</td><td>FWHM Spike</td><td></td><td>84.5</td><td>91.5</td><td>59.5</td><td>60.5</td><td>78.5</td><td>74.9</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td></td><td>85.5</td><td>80.0</td><td>59.5</td><td>55.5</td><td>96.0</td><td>75.3</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td></td><td>77.5</td><td>80.5</td><td>62.5</td><td>80.0</td><td>83.0</td><td>76.7</td></tr><tr><td>SW Delay</td><td>FWHM W/FWHM S</td><td></td><td>54.5</td><td>69.0</td><td>51.5</td><td>68.0</td><td>68.0</td><td>62.2</td></tr><tr><td>Wave Amp</td><td>FWHM W/FWHM D</td><td></td><td>61.0</td><td>58.5</td><td>51.0</td><td>76.0</td><td>68.0</td><td>62.9</td></tr><tr><td>Wave Amp</td><td>FWHM W/FWHM S</td><td></td><td>73.0</td><td>76.5</td><td>51.0</td><td>62.0</td><td>60.0</td><td>64.5</td></tr><tr><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td></td><td>58.0</td><td>76.0</td><td>60.0</td><td>73.5</td><td>62.5</td><td>66.0</td></tr><tr><td>Spike Amp/Wave Amp</td><td>FWHM W/FWHM D</td><td></td><td>52.5</td><td>65.0</td><td>64.5</td><td>79.5</td><td>70.0</td><td>66.3</td></tr><tr><td>SW Delay</td><td>FWHM W/FWHM D</td><td></td><td>62.5</td><td>54.0</td><td>51.5</td><td>81.0</td><td>85.5</td><td>66.9</td></tr><tr><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td></td><td>62.0</td><td>84.0</td><td>63.5</td><td>74.5</td><td>61.5</td><td>69.1</td></tr><tr><td>FWHM Delay</td><td>Spike Amp/Wave Amp</td><td></td><td>65.5</td><td>52.5</td><td>60.5</td><td>80.5</td><td>92.0</td><td>70.2</td></tr><tr><td>Spike Amp</td><td>FWHM W/FWHM D</td><td></td><td>66.0</td><td>84.5</td><td>60.0</td><td>76.5</td><td>68.0</td><td>71.0</td></tr><tr><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td></td><td>70.5</td><td>85.5</td><td>55.5</td><td>59.5</td><td>86.5</td><td>71.5</td></tr><tr><td>Spike Amp</td><td>FWHM W/FWHM S</td><td></td><td>73.5</td><td>89.0</td><td>58.5</td><td>77.0</td><td>61.0</td><td>71.8</td></tr><tr><td>SW Delay</td><td>Spike Amp/Wave Amp</td><td></td><td>64.5</td><td>71.5</td><td>59.5</td><td>82.5</td><td>86.0</td><td>72.8</td></tr><tr><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td></td><td>83.5</td><td>73.5</td><td>55.5</td><td>63.5</td><td>94.0</td><td>74.0</td></tr><tr><td>FWHM Spike</td><td>FWHM W/FWHM D</td><td></td><td>78.5</td><td>88.0</td><td>62.0</td><td>54.5</td><td>87.5</td><td>74.1</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM Wave</td><td>61.0</td><td>68.0</td><td>50.0</td><td>65.5</td><td>69.5</td><td>62.8</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td>FWHM Wave</td><td>75.5</td><td>55.5</td><td>50.5</td><td>64.5</td><td>90.5</td><td>67.3</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM Wave</td><td>53.5</td><td>80.5</td><td>61.5</td><td>80.5</td><td>62.5</td><td>67.7</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>FWHM Wave</td><td>75.0</td><td>61.0</td><td>52.0</td><td>65.0</td><td>86.5</td><td>67.9</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM Wave</td><td>66.5</td><td>82.0</td><td>51.0</td><td>71.5</td><td>70.5</td><td>68.3</td></tr></tbody></table>

(Continued.)
(Continued.)

<table><thead><tr><th>Feature 1</th><th>Feature 2</th><th>Feature 3</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th><th>Acc</th></tr></thead><tbody><tr><td>Wave Amp</td><td>FWHM Spike</td><td>FWHM Wave</td><td>74.5</td><td>80.5</td><td>51.5</td><td>54.0</td><td>86.5</td><td>69.4</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM Delay</td><td>75.0</td><td>58.5</td><td>50.5</td><td>73.0</td><td>91.5</td><td>69.7</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM Delay</td><td>73.5</td><td>58.5</td><td>58.0</td><td>75.5</td><td>91.0</td><td>71.3</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td>FWHM Wave</td><td>78.0</td><td>77.0</td><td>50.5</td><td>65.5</td><td>87.5</td><td>71.7</td></tr><tr><td>Spike Amp</td><td>FWHM Spike</td><td>FWHM Wave</td><td>75.0</td><td>92.5</td><td>56.5</td><td>58.5</td><td>86.5</td><td>73.8</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>FWHM Spike</td><td>85.5</td><td>63.0</td><td>51.5</td><td>74.5</td><td>95.0</td><td>73.9</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td>FWHM Spike</td><td>87.0</td><td>65.0</td><td>50.5</td><td>75.5</td><td>95.5</td><td>74.7</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>SW Delay</td><td>74.5</td><td>74.0</td><td>63.5</td><td>89.5</td><td>77.5</td><td>75.8</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM Delay</td><td>74.5</td><td>72.5</td><td>65.0</td><td>85.5</td><td>82.0</td><td>75.9</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td>FWHM Wave</td><td>87.5</td><td>82.0</td><td>53.5</td><td>64.0</td><td>93.0</td><td>76.0</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM Wave</td><td>78.0</td><td>85.0</td><td>61.0</td><td>75.0</td><td>85.5</td><td>76.9</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM Spike</td><td>88.5</td><td>72.5</td><td>51.5</td><td>81.5</td><td>96.5</td><td>78.1</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM Spike</td><td>83.5</td><td>89.0</td><td>59.5</td><td>83.0</td><td>82.5</td><td>79.5</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM Spike</td><td>86.5</td><td>91.5</td><td>59.5</td><td>79.5</td><td>93.5</td><td>82.1</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM Spike</td><td>91.5</td><td>90.5</td><td>60.0</td><td>78.0</td><td>93.5</td><td>82.7</td></tr><tr><td>SW Delay</td><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td>64.5</td><td>73.5</td><td>54.0</td><td>65.5</td><td>70.5</td><td>65.6</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM W/FWHM D</td><td>60.5</td><td>55.0</td><td>55.0</td><td>81.5</td><td>81.5</td><td>66.7</td></tr><tr><td>Wave Amp</td><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>80.5</td><td>51.5</td><td>51.5</td><td>63.5</td><td>90.0</td><td>67.4</td></tr><tr><td>Wave Amp</td><td>FWHM Spike</td><td>FWHM W/FWHM D</td><td>71.5</td><td>59.0</td><td>51.0</td><td>68.5</td><td>89.0</td><td>67.8</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM W/FWHM D</td><td>63.0</td><td>72.0</td><td>61.0</td><td>83.0</td><td>68.5</td><td>69.5</td></tr><tr><td>SW Delay</td><td>Spike Amp/Wave Amp</td><td>FWHM W/FWHM D</td><td>56.0</td><td>58.5</td><td>61.0</td><td>83.5</td><td>89.5</td><td>69.7</td></tr><tr><td>SW Delay</td><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td>58.5</td><td>79.0</td><td>53.5</td><td>69.5</td><td>89.0</td><td>69.9</td></tr><tr><td>Wave Amp</td><td>SW Delay</td><td>FWHM W/FWHM S</td><td>72.5</td><td>76.0</td><td>50.5</td><td>67.0</td><td>84.0</td><td>70.0</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>75.0</td><td>65.0</td><td>53.0</td><td>71.0</td><td>91.5</td><td>71.1</td></tr><tr><td>SW Delay</td><td>FWHM Delay</td><td>Spike Amp/Wave Amp</td><td>71.5</td><td>55.5</td><td>56.5</td><td>80.5</td><td>91.5</td><td>71.1</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM W/FWHM D</td><td>73.5</td><td>57.0</td><td>63.0</td><td>81.5</td><td>82.0</td><td>71.4</td></tr><tr><td>FWHM Spike</td><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td>74.0</td><td>84.5</td><td>57.5</td><td>53.5</td><td>89.0</td><td>71.7</td></tr><tr><td>FWHM Spike</td><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td>75.5</td><td>87.5</td><td>56.0</td><td>63.0</td><td>79.5</td><td>72.3</td></tr><tr><td>Spike Amp</td><td>Wave Amp</td><td>FWHM W/FWHM S</td><td>77.0</td><td>84.5</td><td>62.5</td><td>79.5</td><td>64.5</td><td>73.6</td></tr><tr><td>Spike Amp</td><td>SW Delay</td><td>FWHM W/FWHM S</td><td>79.0</td><td>85.0</td><td>52.5</td><td>75.5</td><td>86.0</td><td>75.6</td></tr><tr><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td>FWHM W/FWHM D</td><td>71.0</td><td>82.5</td><td>59.0</td><td>79.0</td><td>87.0</td><td>75.7</td></tr><tr><td>FWHM Delay</td><td>FWHM Wave</td><td>Spike Amp/Wave Amp</td><td>74.5</td><td>74.0</td><td>60.0</td><td>81.5</td><td>89.0</td><td>75.8</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td>FWHM W/FWHM D</td><td>89.5</td><td>64.0</td><td>50.0</td><td>81.0</td><td>97.0</td><td>76.3</td></tr><tr><td>Spike Amp</td><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>86.0</td><td>86.0</td><td>59.0</td><td>77.0</td><td>87.5</td><td>79.1</td></tr><tr><td>FWHM Delay</td><td>FWHM W/FWHM S</td><td>Spike Amp/Wave Amp</td><td>80.5</td><td>83.0</td><td>61.0</td><td>80.0</td><td>91.5</td><td>79.2</td></tr><tr><td>Spike Amp</td><td>FWHM Spike</td><td>FWHM W/FWHM D</td><td>89.5</td><td>92.5</td><td>59.5</td><td>76.0</td><td>82.5</td><td>80.0</td></tr><tr><td>FWHM Delay</td><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td>86.5</td><td>85.5</td><td>54.0</td><td>80.5</td><td>95.0</td><td>80.3</td></tr><tr><td>SW Delay</td><td>FWHM Spike</td><td>Spike Amp/Wave Amp</td><td>84.5</td><td>83.5</td><td>57.0</td><td>84.5</td><td>93.5</td><td>80.6</td></tr></tbody></table>

## Appendix D. Model equations and parameters

System equations used for the minimal model reads:

$$ 
\begin{aligned}
\ddot{y}_{\text{PY R}} = & \frac{W_{\text{PY R}}}{\tau_{\text{EPSP}}} S(y_{\text{PY R}}' - C_{\text{PV} \to \text{PYR}} y_{\text{PV}} - C_{\text{SST}_B} \to \text{PYR} y_{\text{SST}_B}) \\
& - C_{\text{SST}_A} \to \text{PYR} y_{\text{SST}_A}) - \frac{2}{\tau_{\text{EPSP}}} \dot{y}_{\text{PY R}} \\
& - \frac{1}{\tau_{\text{EPSP}}^2} y_{\text{PY R}},
\end{aligned}
\quad (4.1a) $$

$$ 
\begin{aligned}
\ddot{y}_{\text{PY R}'} = & \frac{W_{\text{PY R}}}{\tau_{\text{EPSP}}} (p(t) + C_{\text{PY R}'} \to \text{PYR} S(C_{\text{PY R} \to \text{PYR}'} y_{\text{PY R}})) \\
& - \frac{2}{\tau_{\text{EPSP}}} \dot{y}_{\text{PY R}'} - \frac{1}{\tau_{\text{EPSP}}^2} y_{\text{PY R}'},
\end{aligned}
\quad (4.1b) $$

$$ 
\begin{aligned}
\ddot{y}_{\text{PV}} = & \frac{W_{\text{PV}}}{\tau_{\text{IPSP,PV}}} S(C_{\text{PY R} \to \text{PV}} y_{\text{PY R}} - C_{\text{SST}_B} \to \text{PV} y_{\text{SST}_B}) \\
& - \frac{2}{\tau_{\text{IPSP,PV}}} \dot{y}_{\text{PV}} - \frac{1}{\tau_{\text{IPSP,PV}}^2} y_{\text{PV}},
\end{aligned}
\quad (4.1c) $$

$$ 
\begin{aligned}
\ddot{y}_{\text{SST}_B} = & \frac{W_{\text{SST}_B}}{\tau_{\text{IPSP,SST}_B}} S(C_{\text{PY R} \to \text{SST}} y_{\text{SST}_B}) \\
& - \frac{2}{\tau_{\text{IPSP,SST}_B}} \dot{y}_{\text{SST}_B} - \frac{1}{\tau_{\text{IPSP,SST}_B}^2} y_{\text{SST}_B},
\end{aligned}
\quad (4.1d) $$

System variables $y_i$ represent synaptic potentials emitted by PYR, PYR', PV+ (denoted as PV+) and SST+ (denoted by SST$_B$ for the basal contacts and SST$_A$ for the apical contacts) neuronal populations. Parameters $W_i$ and $\tau_i$ correspond to synaptic gain and synaptic time constant, respectively. Coupling coefficients between the subpopulations i and j are denoted by parameters $C_{i \to j}$. Unspecific thalamo-cortical inputs are represented by $p(t) = p_{\text{mean}} + \xi$ with $p_{\text{mean}}$ being the mean and $\xi$ being a random variable following a normal distribution, $N(0, \sigma^2)$. System (4.1a)–(4.1e) is used to simulate the laminar
<table><caption>**Table D1.** Parameter set used to simulate figures 7 and 8. Post-synaptic time are extracted from the literature [67–70].</caption><thead><tr><th>Parameter</th><th>Figure 7</th><th>Figure 8(a)</th><th>Figure 8(b)</th><th>Figure 8(c)</th><th>Figure 8(d)</th><th>Figure 8(e)</th><th>Figure 8(f)</th><th>Figure 8(g)</th><th>Figure 8(h)</th></tr></thead><tbody><tr><td><strong>PSP amplitude</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>W<sub>PY R</sub> (mV)</td><td>8</td><td>8</td><td>8</td><td>8</td><td>8</td><td>8</td><td>8</td><td>8</td><td>16</td></tr><tr><td>W<sub>PV</sub> (mV)</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>W<sub>SST<sub>B</sub></sub> (mV)</td><td>50</td><td>50</td><td>50</td><td>50</td><td>50</td><td>95</td><td>30</td><td>50</td><td>50</td></tr><tr><td>W<sub>SST<sub>A</sub></sub> (mV)</td><td>20</td><td>20</td><td>20</td><td>10</td><td>25</td><td>20</td><td>20</td><td>20</td><td>20</td></tr><tr><td><strong>PSP time constant</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>τ<sub>EPSP</sub> (s)</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.004</td></tr><tr><td>τ<sub>IPSP,PV</sub> (s)</td><td>0.002</td><td>0.002</td><td>0.02</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>τ<sub>IPSP,SST<sub>B</sub></sub> (s)</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.0105</td><td>0.03</td><td>0.02</td><td>0.02</td></tr><tr><td>τ<sub>IPSP,SST<sub>A</sub></sub> (s)</td><td>0.05</td><td>0.05</td><td>0.05</td><td>0.1</td><td>0.04</td><td>0.05</td><td>0.05</td><td>0.05</td><td>0.05</td></tr><tr><td><strong>Coupling coefficient</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sub>PV→PYR</sub></td><td>121</td><td>121</td><td>121</td><td>121</td><td>121</td><td>121</td><td>121</td><td>121</td><td>121</td></tr><tr><td>C<sub>SST<sub>B</sub>→PYR</sub></td><td>26</td><td>50</td><td>20</td><td>26</td><td>26</td><td>26</td><td>26</td><td>26</td><td>26</td></tr><tr><td>C<sub>SST<sub>A</sub>→PYR</sub></td><td>24</td><td>0</td><td>30</td><td>24</td><td>24</td><td>24</td><td>24</td><td>24</td><td>24</td></tr><tr><td>C<sub>PY R→PYR</sub></td><td>108</td><td>108</td><td>108</td><td>108</td><td>108</td><td>108</td><td>108</td><td>108</td><td>108</td></tr><tr><td>C<sub>PY R→PYR'</sub></td><td>135</td><td>135</td><td>135</td><td>135</td><td>135</td><td>135</td><td>135</td><td>135</td><td>135</td></tr><tr><td>C<sub>PY R→PV</sub></td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td></tr><tr><td>C<sub>SST<sub>B</sub>→PV</sub></td><td>40</td><td>40</td><td>40</td><td>40</td><td>40</td><td>40</td><td>40</td><td>40</td><td>40</td></tr><tr><td>C<sub>PY R→SST</sub></td><td>55</td><td>55</td><td>55</td><td>55</td><td>55</td><td>55</td><td>55</td><td>55</td><td>55</td></tr><tr><td><strong>External input</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>p<sub>mean</sub> (Hz)</td><td>90</td><td>90</td><td>90</td><td>90</td><td>90</td><td>90</td><td>90</td><td>90</td><td>90</td></tr><tr><td>σ<sup>2</sup> (Hz<sup>2</sup>)</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td></tr></tbody></table>

<table><caption>**Table D2.** Parameter set used to simulate figure 9. Post-synaptic time are extracted from the literature [67–70].</caption><thead><tr><th>Parameter</th><th>EZ in figure 9</th><th>NEZ in figure 9</th></tr></thead><tbody><tr><td><strong>PSP amplitude</strong></td><td></td><td></td></tr><tr><td>W<sub>PY R</sub> (mV)</td><td>9</td><td>6</td></tr><tr><td>W<sub>PY R''</sub> (mV)</td><td>—</td><td>15</td></tr><tr><td>W<sub>PV</sub> (mV)</td><td>—</td><td>5</td></tr><tr><td>W<sub>SST<sub>B</sub></sub> (mV)</td><td>60</td><td>50</td></tr><tr><td>W<sub>SST<sub>A</sub></sub> (mV)</td><td>35</td><td>20</td></tr><tr><td><strong>PSP time constant</strong></td><td></td><td></td></tr><tr><td>τ<sub>EPSP</sub> (s)</td><td>0.01</td><td>0.01</td></tr><tr><td>τ<sub>EPSP,fast</sub> (s)</td><td>—</td><td>0.01</td></tr><tr><td>τ<sub>EPSP,slow</sub> (s)</td><td>—</td><td>0.05</td></tr><tr><td>τ<sub>IPSP,PV</sub> (s)</td><td>0.002</td><td>0.002</td></tr><tr><td>τ<sub>IPSP,SST<sub>B</sub></sub> (s)</td><td>0.02</td><td>0.02</td></tr><tr><td>τ<sub>IPSP,SST<sub>A</sub></sub> (s)</td><td>0.05</td><td>0.05</td></tr><tr><td><strong>Coupling coefficient</strong></td><td></td><td></td></tr><tr><td>C<sub>PV→PYR</sub></td><td>121</td><td>121</td></tr><tr><td>C<sub>SST<sub>B</sub>→PYR</sub></td><td>27</td><td>16</td></tr><tr><td>C<sub>SST<sub>A</sub>→PYR</sub></td><td>20</td><td>17</td></tr><tr><td>C<sub>PY R→PYR</sub></td><td>108</td><td>100</td></tr><tr><td>C<sub>PY R''→PYR</sub></td><td>—</td><td>2</td></tr><tr><td>C<sub>PY R→PYR'</sub></td><td>135</td><td>135</td></tr><tr><td>C<sub>PY R→PYR''</sub></td><td>—</td><td>100</td></tr><tr><td>C<sub>PY R→PV</sub></td><td>100</td><td>100</td></tr><tr><td>C<sub>PY R''→PV</sub></td><td>—</td><td>0</td></tr><tr><td>C<sub>SST<sub>B</sub>→PV</sub></td><td>40</td><td>40</td></tr><tr><td>C<sub>PY R→SST</sub></td><td>38</td><td>25</td></tr><tr><td>C<sub>PY R''→SST</sub></td><td>—</td><td>3</td></tr><tr><td>C<sub>Ext→PYR</sub></td><td>—</td><td>25</td></tr><tr><td><strong>External input</strong></td><td></td><td></td></tr><tr><td>p<sub>mean</sub> (Hz)</td><td>100</td><td>90</td></tr><tr><td>σ<sup>2</sup> (Hz<sup>2</sup>)</td><td>2</td><td>2</td></tr></tbody></table>
NMM in figures 7 and 8, as well as the EZ model in figure 9. Parameter values are precised in tables D1 and D2.

As for the NEZ model in figure 9 with slow-glutamatergic activation $y_{PY R''}$ and external input $y_{ext}$, system equations read:

$$ 
\begin{aligned}
\ddot{y}_{PY R} = & \frac{W_{PY R}}{\tau_{EPSP,fast}} S(y_{PY R'} + C_{PY R'' \to PYR} y_{PY R''}) \\
& - C_{PV \to PYR} y_{PV} - C_{SST_B \to PYR} y_{SST_B} \\
& - C_{SST_A \to PYR} y_{SST_A} + C_{Ext \to PYR} y_{ext}) \\
& - \frac{2}{\tau_{EPSP,fast}} \dot{y}_{PY R} - \frac{1}{\tau_{EPSP,fast}^2} y_{PY R},
\end{aligned}
\quad (4.2a) $$

$$ 
\begin{aligned}
\ddot{y}_{PY R'} = & \frac{W_{PY R}}{\tau_{EPSP,fast}} (p(t) + C_{PY R' \to PYR} S \\
& \times (C_{PY R \to PYR} y_{PY R}) \\
& - \frac{2}{\tau_{EPSP,fast}} \dot{y}_{PY R'} - \frac{1}{\tau_{EPSP,fast}^2} y_{PY R'}),
\end{aligned}
\quad (4.2b) $$

$$ 
\begin{aligned}
\ddot{y}_{PY R''} = & \frac{W_{PY R''}}{\tau_{EPSP,slow}} S(C_{PY R \to PYR''} y_{PY R''}) \\
& - \frac{2}{\tau_{EPSP,slow}} \dot{y}_{PY R''} - \frac{1}{\tau_{EPSP,slow}^2} y_{PY R''},
\end{aligned}
\quad (4.2c) $$

$$ 
\begin{aligned}
\ddot{y}_{PV} = & \frac{W_{PV}}{\tau_{IPSP,PV}} S(C_{PY R \to PV} y_{PY R} \\
& + C_{PY R'' \to PV} y_{PY R''} - C_{SST_B \to PV} y_{SST_B}) \\
& - \frac{2}{\tau_{IPSP,PV}} \dot{y}_{PV} - \frac{1}{\tau_{IPSP,PV}^2} y_{PV},
\end{aligned}
\quad (4.2d) $$

$$ 
\begin{aligned}
\ddot{y}_{SST_B} = & \frac{W_{SST_B}}{\tau_{IPSP,SST_B}} \times S(C_{PY R \to SST} y_{SST_B} + C_{PY R'' \to SST_B} y_{PY R''}) \\
& - \frac{2}{\tau_{IPSP,SST_B}} \dot{y}_{SST_B} - \frac{1}{\tau_{IPSP,SST_B}^2} y_{SST_B},
\end{aligned}
\quad (4.2e) $$

$$ 
\begin{aligned}
\ddot{y}_{SST_A} = & \frac{W_{SST_A}}{\tau_{IPSP,SST_A}} \times S(C_{PY R \to SST} y_{SST_A} + C_{PY R'' \to SST_A} y_{PY R''}) \\
& - \frac{2}{\tau_{IPSP,SST_A}} \dot{y}_{SST_A} - \frac{1}{\tau_{IPSP,SST_A}^2} y_{SST_A},
\end{aligned}
\quad (4.2f) $$

$$ 
\begin{aligned}
\ddot{y}_{ext} = & \frac{W_{PY R}}{\tau_{EPSP,fast}} (C_{Ext \to PYR} FR_{Ext}) \\
& - \frac{2}{\tau_{EPSP,fast}} \dot{y}_{Ext} - \frac{1}{\tau_{EPSP,fast}^2} y_{Ext}
\end{aligned}
\quad (4.2g) $$

with $FR_{Ext}$ representing the firing rate, hence, the output of the wave-to-rate function $S(.)$ of the pre-synaptic population. Parameter values are given in table D2.

## ORCID iDs

Elif Köksal-Ersöz https://orcid.org/0000-0003-3696-7953

Remo Lazazzera https://orcid.org/0000-0002-4086-3164

Maxime Yochum https://orcid.org/0000-0001-8632-6041

Isabelle Merlet https://orcid.org/0000-0003-2590-6848

Julia Makhalova https://orcid.org/0000-0001-7962-2942

Borja Mercadal https://orcid.org/0000-0002-0941-6796

Roser Sanchez-Todo https://orcid.org/0000-0003-1580-8750

Giulio Ruffini https://orcid.org/0000-0003-3084-0177

Fabrice Bartolomei https://orcid.org/0000-0002-1678-0297

Pascal Benquet https://orcid.org/0000-0003-0217-5886

Fabrice Wendling https://orcid.org/0000-0003-2428-9665

## References

[1] Talairach J and Bancaud J 1966 Lesion, “irritative” zone and epileptogenic focus *Stereotact. Funct. Neurosurg.* 27 91–94

[2] de Curtis M and Avanzini G 2001 Interictal spikes in focal epileptogenesis *Prog. Neurobiol.* 63 541–67

[3] de Curtis M, Jefferys J G R and Avoli M 2012 Interictal epileptiform discharges in partial epilepsy: complex neurobiological mechanisms based on experimental and clinical evidence *Jasper's Basic Mechanisms of the Epilepsies* 4th edn (Bethesda, MD: National Center for Biotechnology Information (US))

[4] Keller C J et al 2010 Heterogeneous neuronal firing patterns during interictal epileptiform discharges in the human cortex *Brain* 133 1668–81

[5] Ulbert I, Heit G, Madsen J, Karmos G and Halgren E 2004 Laminar analysis of human neocortical interictal spike generation and propagation: current source density and multiunit analysis *in vivo Epilepsia* 45 48–56

[6] Hotka M and Kubista H 2019 The paroxysmal depolarization shift in epilepsy research *Int. J. Biochem. Cell Biol.* 107 77–81

[7] Cohen I, Navarro V, Clemenceau S, Baulac M and Miles R 2002 On the origin of interictal activity in human temporal lobe epilepsy *in vitro Science* 298 1418–21

[8] Huberfeld G, de la Prida L M, Pallud J, Cohen I, le van Quyen M, Adam C, Clemenceau S, Baulac M and Miles R 2011 Glutamatergic pre-ictal discharges emerge at the transition to seizure in human epilepsy *Nat. Neurosci.* 14 627–34

[9] Chatrian G E 1974 A glossary of terms most commonly used by clinical electroencephalographers *Electroencephalogr. Clin. Neurophysiol.* 37 538–48

[10] Alarcon G, Garcia Seoane J J, Binnie C D, Martin Miguel M C, Juler J, Polkey C E, Elwes R D and Ortiz Blasco J M 1997 Origin and propagation of interictal discharges in the acute electrocorticogram. Implications for pathophysiology and surgical treatment of temporal lobe epilepsy *Brain* 120 2259–82

[11] Serafini R 2019 Similarities and differences between the interictal epileptiform discharges of green-spikes and red-spikes zones of human neocortex *Clin. Neurophysiol.* 130 396–405
[12] Musolino A, Tournoux P, Missir O and Talairach J 1990 Methodology of “in vivo” anatomical study and stereo-electroencephalographic exploration in brain surgery for epilepsy *J. Neuroradiol.* **17** 67–102 (available at: https://pubmed.ncbi.nlm.nih.gov/2123498/)

[13] Talairach J, Bancaud J, Szikla G, Bonis A, Geier S and Vedrenne C 1974 New approach to the neurosurgery of epilepsy. Stereotaxic methodology and therapeutic results. 1. Introduction and history *Neurochirurgie* **20** 1–240 (available at: https://pubmed.ncbi.nlm.nih.gov/4603712/)

[14] Abd El-Samie F E, Alotaiby T N, Khalid M I, Alshebeili S A and Aldosari S A 2018 A review of EEG and MEG epileptic spike detection algorithms *IEEE Access* **6** 60673–88

[15] Bourien J, Bartolomei F, Bellanger J J, Gavaret M, Chauvel P and Wendling F 2005 A method to identify reproducible subsets of co-activated structures during interictal spikes. Application to intracerebral EEG in temporal lobe epilepsy *Clin. Neurophysiol.* **116** 443–55

[16] Senhadji L, Dillenseger J L, Wendling F, Rocha C and Kinie A 1995 Wavelet analysis of EEG for three-dimensional mapping of epileptic events *Ann. Biomed. Eng.* **23** 543–52

[17] Basseville M and Nikiforov I V 1993 *Detection of Abrupt Changes: Theory and Application* vol 104 (Englewood Cliffs, NJ: Prentice Hall)

[18] Bartolomei F, Chauvel P and Wendling F 2008 Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG *Brain* **131** 1818–30

[19] Chatrian G E 1967 *La Stéréo-électroencéphalographie dans L'épilepsie. Informations Neurophysiopathologiques Apportées par L'investigation Fonctionnelle Stéréotaxique* ed J Bancaud, J Talairach, A Bonis, C Schaub, G Szikla, P Morel and M Bordas-Ferer (Paris: Masson) p 321

[20] Allen P J, Fish D R and Smith S J M 1992 Very high-frequency rhythmic activity during SEEG suppression in frontal lobe epilepsy *Electroencephalogr. Clin. Neurophysiol.* **82** 155–9

[21] Alarcon G, Binnie C D, Elwes R D C and Polkey C E 1995 Power spectrum and intracranial EEG patterns at seizure onset in partial epilepsy *Electroencephalogr. Clin. Neurophysiol.* **94** 326–37

[22] Wendling F, Bartolomei F, Bellanger J J, Bourien J and Chauvel P 2003 Epileptic fast intracerebral EEG activity: evidence for spatial decorrelation at seizure onset *Brain* **126** 1449–59

[23] Pedregosa F et al 2011 Scikit-learn: machine learning in Python *J. Mach. Learn. Res.* **12** 2825–30

[24] da Silva F L, Hoeks A, Smits H and Zetterberg L H 1974 Model of brain rhythmic activity: the alpha-rhythm of the thalamus *Kybernetik* **15** 27–37

[25] Freeman W J 1987 Simulation of chaotic EEG patterns with a dynamic model of the olfactory system *Biol. Cybern.* **56** 139–50

[26] Jansen B H and Rit V G 1995 Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns *Biol. Cybern.* **73** 357–66

[27] Wendling F, Bartolomei F, Bellanger J J and Chauvel P 2002 Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition *Eur. J. Neurosci.* **15** 1499–508

[28] Köksal-Ersöz E, Modolo J, Bartolomei F and Wendling F 2020 Neural mass modeling of slow-fast dynamics of seizure initiation and abortion *PLoS Comput. Biol.* **16** e1008430

[29] Bhattacharya B S, Coyle D and Maguire L P 2011 A thalamo-cortico-thalamic neural mass model to study alpha rhythms in Alzheimer's disease *Neural Netw.* **24** 631–45

[30] Ruffini G, Sanchez-Todo R, Dubreuil L, Salvador R, Pinotsis D, Miller E K, Wendling F, Santarnecchi E and Bastos A 2020 P118 a biophysically realistic laminar neural mass modeling framework for transcranial current stimulation *Clin. Neurophysiol.* **131** e78–e97

[31] Sanchez-Todo R, Bastos A, Lopez-Sola E, Mercadal B, Santarnecchi E, Miller E K, Deco G and Ruffini G 2022 A

physical neural mass model framework for the analysis of oscillatory generators from laminar electrophysiological recordings *bioRxiv Preprint* (https://doi.org/10.1101/ 2022.07.19.500618) (Posted online 27 August 2022)

[32] Lopez-Sola E et al 2021 A personalizable autonomous neural mass model of epileptic seizures *J. Neural Eng.* **19** 055002

[33] DeFelipe J, Alonso-Nanclares L and Arellano J I 2002 Microstructure of the neocortex: comparative aspects *J. Neurocytol.* **31** 299–316

[34] Cosandier-Rimélé D, Badier J M, Chauvel P and Wendling F 2007 A physiologically plausible spatio-temporal model for EEG signals recorded with intracerebral electrodes in human partial epilepsy *IEEE Trans. Biomed. Eng.* **54** 380–88

[35] Buzsáki G, Anastassiou C A and Koch C 2012 The origin of extracellular fields and currents—EEG, ECoG, LFP and spikes *Nat. Rev. Neurosci.* **13** 407–20

[36] Speckmann E J and Elger C E 1999 Introduction to the neurophysiological basis of the EEG and DC potentials *Electroencephalography: Basic Principals, Clinical Applications and Related Fields* 4th edn (Baltimore: Williams and Wilkins)

[37] Wang Y, Toledo-Rodriguez M, Gupta A, Wu C, Silberberg G, Luo J and Markram H 2004 Anatomical, physiological and molecular properties of Martinotti cells in the somatosensory cortex of the juvenile rat *J. Physiol.* **561** 65–90

[38] Markram H et al 2015 Reconstruction and simulation of neocortical microcircuitry *Cell* **163** 456–92

[39] Eyal G, Verhoog M B, Testa-Silva G, Deitcher Y, Benavides-Piccione R, DeFelipe J, De Kock C P J, Mansvelder H D and Segev I 2018 Human cortical pyramidal neurons: from spines to spikes via models *Front. Cell Neurosci.* **12** 181

[40] Dinocourt C, Petanjek Z, Freund T F, Ben-Ari Y and Esclapez M 2003 Loss of interneurons innervating pyramidal cell dendrites and axon initial segments in the CA1 region of the hippocampus following pilocarpine-induced seizures *J. Comp. Neurol.* **459** 407–25

[41] Drexel M, Rahimi S and Sperk G 2022 Silencing of hippocampal somatostatin interneurons induces recurrent spontaneous limbic seizures in mice *Neuroscience* **487** 155–65

[42] Huberfeld G, Wittner L, Clemenceau S, Baulac M, Kaila K, Miles R and Rivera C 2007 Perturbed chloride homeostasis and GABAergic signaling in human temporal lobe epilepsy *J. Neurosci.* **27** 9866–73

[43] Burman R J, Rosch R E, Wilmhurst J M, Sen A, Ramantani G, Akerman C J and Raimondo J V 2022 Why won't it stop? The dynamics of benzodiazepine resistance in status epilepticus *Nat. Rev. Neurol.* **18** 428–41

[44] Abegg M H, Savic N, Ehrengruber M U, McKinney R A and Gåhwiler B H 2004 Epileptiform activity in rat hippocampus strengthens excitatory synapses *J. Physiol.* **554** 439–48

[45] Naylor D E, Liu H, Niquet J and Wasterlain C G 2013 Rapid surface accumulation of NMDA receptors increases glutamatergic excitation during status epilepticus *Neurobiol. Dis.* **54** 225–38

[46] Hardie J B and Pearce R A 2006 Active and passive membrane properties and intrinsic kinetics shape synaptic inhibition in hippocampal CA1 pyramidal neurons *J. Neurosci.* **26** 8559–69

[47] Czeiger D and White E L 1993 Synapses of extrinsic and intrinsic origin made by callosal projection neurons in mouse visual cortex *J. Comp. Neurol.* **330** 502–13

[48] Thomson A M, West D C, Hahn J and Deuchars J 1996 Single axon IPSPs elicited in pyramidal cells by three classes of interneurones in slices of rat neocortex *J. Physiol.* **496** 81–102

[49] Deisz R A 1999 The GABA<sub>B</sub> receptor antagonist CGP 55845A reduces presynaptic GABA<sub>B</sub> actions in neocortical neurons of the rat *in vitro Neuroscience* **93** 1241–9

[50] Nettleton J S and Spain W J 2000 Linear to supralinear summation of AMPA-mediated EPSPs in neocortical pyramidal neurons *J. Neurophysiol.* **83** 3310–22
[51] Côme E, Marques X, Poncer J C and Lévi S 2020 KCC2 membrane diffusion tunes neuronal chloride homeostasis *Neuropharmacology* **169** 107571

[52] Serafini R and Loeb J A 2015 Enhanced slow waves at the periphery of human epileptic foci *Clin. Neurophysiol.* **126** 1117–23

[53] Köksal-Ersöz E and Wendling F 2021 Canard solutions in neural mass models: consequences on critical regimes *J. Math. Neurosci.* **11** 11

[54] Tremblay R, Lee S and Rudy B 2016 GABAergic interneurons in the neocortex: from cellular properties to circuits *Neuron* **91** 260–92

[55] Urban-Ciecko J, Fanselow E E and Barth A L 2015 Neocortical somatostatin neurons reversibly silence excitatory transmission via GABAb receptors *Curr. Biol.* **25** 722–31

[56] Obermayer J, Heistek T S, Kerkhofs A, Goriounova N A, Kroon T, Baayen J C, Idema S, Testa-Silva G, Couey J J and Mansvelder H D 2018 Lateral inhibition by Martinotti interneurons is facilitated by cholinergic inputs in human and mouse neocortex *Nat. Commun.* **9** 4101

[57] Hagger-Vaughan N and Storm J F 2019 Synergy of glutamatergic and cholinergic modulation induces plateau potentials in hippocampal OLM interneurons *Front. Cell Neurosci.* **13** 508

[58] Kroon T et al 2019 Group I mGluR-mediated activation of Martinotti cells inhibits local cortical circuitry in human cortex *Front. Cell Neurosci.* **13** 315

[59] Perreault P and Avoli M 1992 4-aminopyridine-induced epileptiform activity and a GABA-mediated long-lasting depolarization in the rat hippocampus *J. Neurosci.* **12** 104–15

[60] Avoli M, Louvel J, Drapeau C, Pumain R and Kurcewicz I 1995 GABAA-mediated inhibition and *in vitro* epileptogenesis in the human neocortex *J. Neurophysiol.* **73** 468–84

[61] Avoli M, Barbarosie M, Lücke A, Nagao T, Lopantsev V and Köhling R 1996 Synchronous GABA-mediated potentials and epileptiform discharges in the rat limbic system *in vitro* *J. Neurosci.* **16** 3912–24

[62] Mercadal B, Lopez-Sola E, Galan-Gadea A, Al-Harrach M, Sanchez-Todo R, Salvador R, Bartolomei F, Wendling F and Ruffini G 2022 Towards a mesoscale physical modeling framework for stereotactic-EEG recordings *bioRxiv Preprint* (https://doi.org/10.1101/2022.07.06.498826) (Posted online 08 July 2022)

[63] Aeed F, Shnitzer T, Talmon R and Schiller Y 2020 Layer- and cell-specific recruitment dynamics during epileptic seizures *in vivo Ann. Neurol.* **87** 97–115

[64] Hall S, Hunt M, Simon A, Cunningham L G, Carracedo L M, Schofield I S, Forsyth R, Traub R D and Whittington M A 2015 Unbalanced peptidergic inhibition in superficial neocortex underlies spike and wave seizure activity *J. Neurosci.* **35** 9302–14

[65] Hall S P, Traub R D, Adams N E, Cunningham M O, Schofield I, Jenkins A J and Whittington M A 2018 Enhanced interlaminar excitation or reduced superficial layer inhibition in neocortex generates different spike-and-wave-like electrographic events *in vitro* *J. Neurophysiol.* **119** 49–61

[66] Goodfellow M 2011 Spatio-temporal modelling and analysis of epileptiform EEG *PhD Thesis* The University of Manchester

[67] Bacci A, Rudolph U, Huguenard J R and Prince D A 2003 Major differences in inhibitory synaptic transmission onto two neocortical interneuron subclasses *J. Neurosci.* **23** 9664–74

[68] Deleuze C, Bhumbra G S, Pazienti A, Lourenço J, Mailhes C, Aguirre A, Beato M and Bacci A 2019 Strong preference for autaptic self-connectivity of neocortical PV interneurons facilitates their tuning to γ-oscillations *PLoS Biol.* **17** e3000419

[69] Karnani M M, Jackson J, Ayzenshtat I, Tucciarone J, Manoocheri K, Snider W G and Yuste R 2016 Cooperative subnetworks of molecularly similar interneurons in mouse neocortex *Neuron* **90** 86–100

[70] Shamas M, Benquet P, Merlet I, Khalil M, El Falou W, Nica A and Wendling F 2018 On the origin of epileptic high frequency oscillations observed on clinical electrodes *Clin. Neurophysiol.* **129** 829–41