Whole-brain simulation of interictal epileptic discharges for patient-specific interpretation of interictal SEEG data

In patients with refractory epilepsy, the clinical interpretation of stereoelectroencephalographic (SEEG) signals is crucial to delineate the epileptogenic network that should be targeted by surgery. We propose a pipeline of patient-specific computational modeling of interictal epileptic activity to improve the definition of regions of interest. Comparison between the computationally defined regions of interest and the resected region confirmed the efficiency of the pipeline. This result suggests that computational modeling can be used to reconstruct signals and aid clinical interpretation.


Introduction
Stereoelectroencephalographic (SEEG) signals recorded in patients with drug-resistant focal epilepsy show distinct activity patterns during interictal, preictal, and ictal periods.These patterns are analyzed to define the epileptogenic zone and propagation zone networks [3].In the clinical routine, this analysis is usually performed visually.Although recorded within epileptic tissue, SEEG signals reflect the projection of epileptic activity from multiple sources onto the multiple contacts of the intracerebral electrode, which complicates their visual interpretation.Previous studies have suggested that computational modeling of source activity can improve the source localization problem of transcranial recordings [6,7].This short communication is a proof of concept illustrating that patient specific modeling of head anatomy and source activity with a neurophysiologically plausible model combined with source localization could improve the interpretation of SEEG recorded interictal epileptic discharges (IEDs) and provide additional information to the clinical interpretation of SEEG signals.Results suggest that patient specific modeling can refine the definition of the Region Of Interest (ROI) responsible for IEDs.

Patient data
The patient was a 14-year-old right-handed girl diagnosed with left occipital, MRI (magnetic resonance imaging)-negative focal cortical dysplasia (FCD).The patient experienced weekly seizures comprising isolated auras characterized by simple visual hallucinations, as well as focal seizures with alteration of awareness (FAES) starting with the same subjective semiology and frequently evolving to bilateral tonic-clonic seizures.
The epileptogenic zone network was defined by expert clinicians (FB, JM) using both visual and quantitative SEEG-signal analysis of interictal and ictal activity based on the epileptogenicity index method [1].The epileptogenic zone was located within the left lateral and polar occipital cortex (contacts GL '4-8) where the ictal and interictal epileptic activity was recorded.Surgery was restricted to the resection of the lateral part of the epileptogenic zone (contacts GL'6-8, including the cortical dysplastic tissue) due to the functional constraints.It led to significant post-surgical improvement in seizures (rare disabling seizures, Engel class II).At the last follow-up 12 years post-surgery, some rare auras and FIAS persisted, with more than 90 % reduction in seizure frequency compared to baseline, whereas no bilateral tonic-clonic seizures occurred after surgery.Histological examination revealed the presence of an FCD type I in the whole resection specimen.

Simulation of SEEG signals on the whole brain model
The patient-specific modeling pipeline (Fig. 2) includes construction of the patient-specific head model and simulation of interictal epileptic activity at the level of SEEG contacts.
For the construction of the patient-specific head model, the DICOM data obtained from MRI scans were converted into NIfTI (Neuroimaging Informatics Technology Initiative) format and then segmented using FreeSurfer [8].The resulting cortical mesh was imported to Brainstorm [21], and down sampled to 15,000 vertices.The neocortex was then subdivided into 144 neocortical regions according to the Virtual Epileptic Patient (VEP) atlas [22].Streamline tractography specifically from the patient's diffusion-weighted MRI (DW-MRI) was used to estimate structural connectivity (number of streamlines) and delay (length of streamlines) matrices [5,19].
IEDs were detected visually from the patient's SEEG data.Detected events were averaged at their peaks.To localize the source of averaged intracerebral IEDs we first calculated a leadfield matrix A which represents the contribution of dipole at each vertex of the cortical mesh (15,000) at the level of each 117 SEEG electrode contacts.The coordinates of the central point of SEEG contacts were estimated with the GARDEL tool [18] from co-registration of the patient's MRI on her CT scan acquired with electrodes in place.The matrix A (117 × 15,000) was then generated from a 1-layer (inner skull, conductivity value of 0.33 S/m) Boundary Element Method (BEM) model using OpenMEEG [4,11] implemented in Brainstorm software [21].Since the weighted minimum-norm estimation (wMNE) is adapted to localize cortical activity of small size, it was used for obtaining the first approximation of ROIs from the peak of the averaged IEDs.We used the default parameters of the unconstrained dipole orientations with depth weight order 0.5, maximal amount 10, noise covariance regularization 0.1 and signal-to-noise ratio 3, output mode as inverse kernel.The estimated ROIs were inserted into the brain atlas by subtracting them from the regions of the VEP atlas.
Forward modeling corresponds to obtaining a simulated SEEG signal from the simulated source activity.The source activity was simulated by using a layered neural mass model (NMM) of the neocortex [13], which was schematized in Fig. 2 (panel "Forward Problem").This model formulates the average activity of the excitatory and inhibitory cell types of the neocortex and their synaptic interactions.It assumes that only the synaptic inputs on the layer V pyramidal cell subpopulation contribute to SEEG signals.Two synaptic positions on the layer V pyramidal cell subpopulation were considered, i.e. layer I for apical projections, and layer V for basal dendritic and perisomatic projections.The extracellular potentials generated by the postsynaptic currents at apical and basal regions were projected onto an artificial electrode contact which is placed 1 mm away from the source activity by using the two-monopole approximation [15].This mesoscale source activity is then considered during the solution to the forward problem (see below).For details of the model architecture, reader may refer to [13,24].At macroscale level each region of the brain atlas was associated with a layered NMM interconnected via the structural connectivity and delay matrices.The parameters of the NMM were identified to fit the averaged IED morphology observed in the real SEEG recordings, not the temporal estimation of the wMNE or the background activity.Model equations and parameters are given in Supplementary Material.
There are two ways of triggering an IED in the computational model.The first one is applying a Gaussian input with a large variance.However, in this case it is not possible to control for either temporal difference between the events in the network or their morphologies.The second approach is perturbing the system with a brief pulse signal, by which we can control the temporal difference between the events.In this study, we opted for the second method to control the temporality and to have less influence on the morphology of the IEDs.
From this model, the same source activity was assigned to all vertices within the brain region associated with an NMM for obtaining the extended source activity.An active (i.e.epileptic) ROI contained dipoles assigned with the source activity signals that simulate an IED.The dipoles of an inactive (i.e.non-epileptic) region were assigned with source activity signals simulating background activity.At the end of this procedure, the simulated source activity defined the extended source matrix S of 15,000 x T, with T being the size of the time vector of the simulated activity.The X (117 x T) matrix of simulated SEEG signals was then obtained by X = A.S.
The patch optimization was performed manually.The ROI provided by the wNME method was considered as the initial guess of the patch size and its location.The patch size a priori determines the amplitude of the reconstructed signal.The amplitude also depends on the patch location: a reconstructed signal from a patch of size N located inside a sulcus would have a smaller amplitude than a reconstructed signal from a patch of size N located on a gyrus.Therefore, we worked on the size and position of the patch to refine the forward solution and we modified

Results
Two distinct types of IEDs were detected visually on different SEEG contacts.The first type of IED was a spike involving the GL'5-GL'6 contacts exploring left occipital cortex, and was preceded by a brief spike 10 ms earlier with the same polarity in the GL'4-GL'5 contacts (Fig. 3A).A polarity inversion was observed between the contacts GL'5-GL'6 and GL'6-GL'7.The spike activity was localized from an average of 20 IEDs by the wMNE in the left lateral occipital and in the left lingual sulcus (basal occipital cortex).The SEEG simulated from the source activities assigned to the wMNE-defined ROIs showed the largest component on GL'4-GL'5 but with a polarity inversion on adjacent contacts (Fig. 3B).The wMNE-defined ROIs were re-delineated by downsizing the one in the left sulcus and by downsizing the one in the left lateral occipital but adding regions which were not included by the wMNE.The new simulated signal showed the same polarity inversion as in the real SEEG, with a 10 ms difference between the events (Fig. 3C).
The second type of IED was a spike-wave with maximal spike amplitude in GL'2-GL'3 and a discrete polarity inversion on adjacent contacts during the first peak of the spike-wave.The transient activity on GL'1-GL'2 preceded the spike-and-wave by 15 ms.The wMNE solution obtained from an average of 20 IEDs was located in the left medial occipital region (lingual gyrus).The SEEG signal simulated by assigning the source activity to the wMNE-defined ROI yielded a maximal spike amplitude on GL'1-GL'2 (Fig. 3E).We optimized the ROI by enlarging the wMNE-defined ROI and dividing it into two (one small ROI is surrounded by a larger ROI).We associated distinct source activities with each region.The new simulated signal showed the maximal spike amplitude on GL'2-GL'3, followed the correct polarity inversions on adjacent contacts and the time difference (Fig. 3F).
The optimization of ROI for forward modeling of IEDs above followed different delineation approaches.The wMNE-defined region was downsized for modeling the IED of spike morphology, whereas it was divided for modeling the IEDs of spike-wave morphology.We exemplified the effect of the ROI definition on the forward solution for these different delineation approaches in Supplementary Figures 1 and 2. The resulting signals do not match the inversion of the polarity of the SEEG recordings.While it is hard to associate one approach with an IED type, this investigation shows again the effect of ROI definition on the forward solution.
The surgically resected zone corresponding to the left lateral occipital (Fig. 3B-F) covered only one of the three regions defined by the model.It is worth mentioning that the surgical resection targeted the clinically defined epileptogenic zone network, which did not include the lingual gyrus and the lingual sulcus.Despite the presence of IEDs, these medial occipital regions were considered by clinicians as a part of the propagation zone.

Discussion
We have presented a pipeline to define cortical ROIs responsible for IED generation by combining SEEG source localization, realistic source modeling, and SEEG forward solution in a patient-specific anatomical model.The source localization from intracerebral recordings is refined by using a computational model that can finely simulate IEDs respecting the polarity inversions on SEEG electrodes, which are determinant for the localization of epileptogenic regions.
In the present case, there was also an incomplete concordance between the clinically defined epileptogenic regions and the IEDgenerating regions simulated by the model.This is not surprising given that the concordance between the regions showing maximal interictal spiking activity and those with maximal epileptogenicity is known to be incomplete in about 25 % of patients with FCD [2].Nonetheless, one of three IED-generating regions identified by the model is in agreement with the resected epileptogenic zone.The other two regions spared by the resection might correspond to the remaining epileptogenic regions responsible for the persisting seizures.
The clinical interpretation of the epileptogenic zone is based on ictal analysis.It has been previously shown that ictal patterns can be simulated using simpler NMM with similar formulations [10,16,23].Our recent study models ictal dynamics in small networks [12].Therefore, the computational approach introduced here can be used to estimate ROIs from SEEG ictal data.
Our method of ROI localization is based on matching averaged IED morphology and polarity inversion of recorded and simulated signals, for which the patch location and its size were also crucial.Modeling IED morphologies allows to decipher the pathophysiological mechanisms mediating IEDs, which can be of clinical relevance [9] and adds information to ictal analysis.Our whole-brain modeling approach could be considered for studying the acute and long-term effects of different stimulation techniques, such as direct cortical stimulation and radiofrequency thermocoagulation.The manual estimation of the model parameters and patches is currently a limitation, which can be overcome by the development of an automatic model fitting pipeline, especially for modeling individual events which would allow the study of study propagation processes during the course of the interictal spike.The overall framework can be adapted to train artificial neural networks for improving SEEG source imaging, as suggested for transcranial imaging techniques [20].
Among different brain parcellation methods, the VEP atlas [22] was used for this methodological study since it has been used in clinical practice at the Epilepsy Unit of the La Timone Hospital in Marseille where the patient was hospitalized, and it was used to determine the patient specific connectivity matrices.That being said, the choice of specific parcellation method is less important for our objective since it is used to refine the definition of the regions of interest (epileptic regions) on the cortical surface.In this patient, the interictal epileptic activity was local and polarity inversion was intact between individual events.Therefore, modeling the epileptic regions for averaged IEDs was sufficient for the forward solution.Nevertheless, we detected and studied two types of IEDs in this patient, which is as important as studying the temporal and spatial differences of one type of event.The non-active (non-epileptic) regions were assigned with a background activity.
We considered a BEM based three-layer homogeneous isotropic head model, although a finite element method (FEM) is more appropriate for SEEG-based analysis.In the future, results could be improved if a FEM is used for simulating SEEG signals.However, it would require using volumes of interest instead of ROIs, as well as a volumetric parcellation of the brain atlas, which is not the case for the surface-based VEP atlas.Finally, the FCD subtypes present specific and distinct alternations of neuronal populations, as well as the architectural and cytoarchitectural abnormalities of cortical layers [14,17].Such pathophysiological abnormalities are difficult to model and might be the current limitation for proposing a more specific modeling approach integrating the underlying etiology.
To conclude, this study suggests that a combination of anatomical and computational models can be successfully used to guide the localization of epileptogenic regions.Although further validation in a large patient cohort is mandatory, this modeling approach is particularly relevant for designing surgical interventions and neuromodulation protocols.Fig. 3. Two distinct types of IEDs were detected on the GL' electrode.(A) Type-I IED (right panel) shows a biphasic spike morphology.The wMNE solution corresponds to the peak of the averaged negative phase of IEDs which is localized in the left lateral occipital and left lingual sulcus.(B) Simulated source activities for the active ROIs (ROI-1 and ROI-2 signals on the left panel) were assigned for the regions defined by the wMNE method, i.e.ROI-1 for the left lateral occipital (blue region) and ROI-2 for the left lingual sulcus (orange region) (middle panels).The background (BKG) source activity was assigned to the inactive regions.The simulated SEEG showed the largest component on GL'4-GL'5 with an inversion of polarity on GL'5-GL'6 and a 10 ms difference between the signals on GL'5-GL'6 and GL'6-GL'7 (right panel).(C) The same source activities for the active and inactive regions were associated with the manually optimized ROIs (middle panels).The new simulated SEEG showed the same polarity on GL'4-GL'5 and GL'5-GL'6 with a 10 ms difference between the peaks and a polarity inversion between GL'5-GL'6 and GL'6-GL'7 (right panel).(D) Type-II IED (right panel) shows a spike-wave morphology.The wMNE solution corresponds to the peak of the averaged positive peak of IEDs which is localized in the left lingual gyrus.(E) Simulated source activities for the active ROI (ROI signal on the left panel) were assigned for the ROI (green patch on the middle panel) defined by the wMNE method.The background (BKG) source activity was assigned to the inactive regions.The simulated SEEG showed the largest component on GL'1-GL'2 with an inversion of polarity on GL'2-GL'3.No activity appeared on GL'3-GL'4 contacts (right panel).(F) The wMNE-derived ROI was optimized by defining a small ROI (ROI-1, yellow patch) and a larger ROI (ROI-2, green patch) (middle panels).A source signal of spike-and-wave morphology was associated with ROI-2 (green ROI-2 signal on the left panel).A source signal that appeared 15 ms before the main spike-and-wave was associated with ROI-1 (yellow ROI-1 signal on the left panel).The new simulated SEEG followed the inverted polarities on GL'1-GL'2, GL'2-GL'3, and GL'3-GL'4, as well as the 15 ms time difference during the initial phase of the SEEG signals on GL'1-GL'2 and GL'2-GL'3 (right panel).Red zones delineated in (B), (C), (E), and (F) show the resected zone in the surgery.Insets zoom into the active regions of interest.

Fig. 1 .
Fig. 1.Posterior view (A) and lateral view (B) of the cortex with SEEG electrodes.

Fig. 2 .
Fig. 2. Pipeline for patient specific intracerebral IED simulation.Patient data including 3D-T1 MRI, DW MRI, CT-scan and SEEG recordings are processed for patient specific head modeling, connectivity matrix, IED detection and categorization.Following the solution to the inverse problem, ROIs were defined and inserted into the atlas.SEEG signals were simulated using a neocortical neural mass model.An iterative process of redefining ROIs and improving the SEEG morphology was performed to reduce the difference between the real and simulated SEEG signals.For the sake of clarity, only CU', FGA' and GL' SEEG electrodes are shown.