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

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.


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 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.
propagation. This categorization is used for deciding the best therapeutic strategy.
Numerous studies have been dedicated to interictal paroxysmal events in humans [1]. Interictal epileptiform discharges (IEDs) observed in SEEG recordings demonstrate a wide range of morphology, ranging from simple monophasic epileptic spikes to more complex multiphasic transient events. Generation of IEDs in partial epilepsies is commonly ascribed to enhanced excitatory interactions within glutamatergic neural networks. However, both human data and animal models have supported the view that gamma-aminobutyric acidergic (GABAergic) signaling does play a significant role [2,3]. The presence of heterogeneous neural firing patterns during IEDs in epileptic and non-epileptic zones suggests a multitude of agents taking part in IED generation [4], including paroxysmal depolarization shifts (PDSs) [5,6] and depolarizing GABAergic actions [7,8]. So far, mechanisms leading to the generation and propagation of SEEG-recorded IEDs remain unclear, and their clinical value for predicting the epileptogenicity 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 spikewaves (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 epilepsy (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 noninvasive methods: e.g. electroencephalogram (EEG), MRI, and computed tomography scan. Figure 1 displays 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 classified according to their morphologies, durations, and generation sites. Their role in interictal epileptogenic 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 laminar neural mass models (NMMs).

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.

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-tonoise 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.

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.

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][21][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).  Patient  P1  P2  P3  P4  P5   Ictal records  2  4  1  4  1  Average ictal records duration (min)  25  38  60  22  30  Inter-ictal records  5  2  9  6  5  Average inter-ictal records duration (min)  52  60  56  53  60  SEEG contacts  93  116  95  116  62  Detected SWs per bipolar channel over all  interictal records   3850  555  7103  3987  4269   Total detected SWs on interictal records  358 024  64 402  674 809  462 518  264 683   Table 2. SEEG experts selected one bipolar channel located in the EZ and one positioned in the NEZ for each patient.

SWs mean waveform computation
The last step consisted of computing the mean waveforms 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 zscore. 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.

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), 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. 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 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.

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 kmeans 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.

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][28][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: A second-order differential equation relates the average presynaptic firing rate to post-synaptic potential (PSP):ÿ where W represents the average synaptic gain, and τ w is the synaptic time constant lumping the rise and decay times (assumed equal), and axonal delays. The time constant parameter τ 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 σ = 0.3 × 10 −3 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 ∞ = 0 is then the difference between the source and the sink, 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), where η = 10 −3 S is a conversion factor relating the PSP to post-synaptic current [31,32]. In reality, η 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 s and the same noise vector following a normal distribution.

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

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.

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 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  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.

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 kmeans 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 interpatient differences in the binary classification task.  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.

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.

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][31][32].  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   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 glutamatergic 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.  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][47][48][49][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.

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 5 and enhanced glutamatergic 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 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 longlasting 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.

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.

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 extracted 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 correctness of extracted features but also limited the number of samples in the set. Among the extracted features, 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 waveforms and the quantitative results of unsupervised classification (tables in appendices B and C), the following 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 correlated 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 features, we have reached an accuracy of 72.6% in separating 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 accuracy 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 combination of primary features, and their use was exploited to discover if they could lead to better classification performances. Finally, these features have not represented a real advantage for classification.
Some comments about patients' specific features distribution (figure in appendix A) are worth mentioning. Particularly for patients P2 and P3, it is visually difficult to affirm that the EZ and NEZ distributions of some features are different. We can consider two causes: the smaller number of data points (100 points per distribution) and the patient-specific morphology. In particular, patient P2 feature distributions 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 activity. Classification has shown poor discrimination performances 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 comparison with literature, Serafini and Loeb [52] have shown that on human SEEG, IEDs at the center of epileptic 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 distribution 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 morphology 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 discriminating features could be present in other different morphologies worth to be analyzed. The presented results are consistent for patients with a specific disorder (TLE). A larger cohort of patients would consolidate 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 directionality of information flow among brain areas.

Spatially distributed synaptic interaction for the SW morphology
In order to model the IEDs measured by SEEG we have used a neurophysiologically-plausible laminar NMM. The synaptic interaction at different layers with synaptic time constants accounting for local dendritic delays has permitted to simulate realistic SW morphologies. While the simple model was responding to the SW morphology, slow glutamatergic signaling was included to account for the temporal 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 epileptogenicity 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 A -mediated inhibition but also slow GABA Bmediated 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 nonsynaptic 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][60][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.

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-towave 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.

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.

Appendix D. Model equations and parameters
System equations used for the minimal model reads: Parameters W i and τ i correspond to synaptic gain and synaptic time constant, respectively. Coupling coefficients between the subpopulations i and j are denoted by parameters C i→j . Unspecific thalamocortical inputs are represented by p(t) = p mean + ξ with p mean being the mean and ξ being a random variable following a normal distribution, N(0, σ 2 ). System (4.1a)-(4.1e) is used to simulate the laminar   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 slowglutamatergic activation y PY R ′ ′ and external input y ext , system equations read: