Machine learning approach to detect focal-onset seizures in the human anterior nucleus of the thalamus

Seizures (N = 79 from 10 patients) and interictal SEEG data (N = 158 segments, i.e. 79 × 2 states from sleep and awake) were obtained from patients with drug-resistant temporal lobe epilepsy (TLE) who have undergone SEEG exploration for epilepsy surgery. The demographic details of the study participants are detailed in (table 1). The Institutional Review Board of the University of Alabama at Birmingham approved the study for recording LFPs from the thalamus during SEEG exploration (IRB-170323005). Before the surgery, informed consent detailing the thalamic implantation for research was obtained from the participants. The multi-step consenting process has been described in detail in a previous study [28]. To mitigate the risk of increased complications from an additional depth electrode placement for research purposes only, the surgeon modified the trajectory of a clinically indicated multi-contact depth electrode to sample LFPs from the insular operculum region (superficial contacts) and the ANT (deeper contacts) (figure 1(A)). Only one ANT was sampled per participant. The rationale for selecting the ANT was based on its—a) structural connectivity with the limbic network [21]; b) preclinical and clinical studies demonstrating recruitment of ANT in TLE [18, 29]; and c) the potentials of translating the knowledge gained in the study towards the development of a closed-loop DBS therapy. Only patients with confirmed electrode localization in the ANT were included in the study

Zoom In Zoom Out Reset image size

Robotic assistance (ROSA device, Medtech, Syracuse, NY) was used to plan optimal and safe trajectories for SEEG multi-contact electrode implantation (12–16 contacts per depth electrode, 2 mm contact length, 0.8 mm contact diameter, 1.5 mm inter-contact distance, PMT^® Corporation, Chanhassen, MN). Intracranial video-EEG was recorded using Natus Quantum (Natus Medical Incorporated, Pleasanton, CA, sampling rate 2048 Hz) (figure 1). Signals were referenced to a common extracranial electrode placed posteriorly in the occiput near the hairline.

The pipeline and the accuracy of targeting thalamic subnuclei during SEEG have been reported in great detail previously by our group [28]. Electrode contacts in the ANT were localized using Lead-DBS v2 software [30] and the cortical contacts using iElectrodes [31]. Preoperative MRI and postoperative CT scans were linearly co-registered using Advanced Normalization Tools followed by refinement with brain shift correction to improve the registration of subcortical structures Jenkinson [32]. Both the images were normalized to ICBM 2009b NLIN asymmetric space using the symmetric diffeomorphic image registration. The cortical regions of all the implanted SEEG contacts were confirmed with AAL2 atlas [33], and the thalamic nuclei were identified using the mean histological thalamic Morel's atlas (figures 1(A) and (B)) [34]. The trajectory of the ANT implantation passed through the precentral gyrus, pars opercularis frontalis, insula, and putamen before ending in the ANT (figure 1(A)) [28].

Seizure onset and offset in the cortex were annotated by a board-certified epileptologist (SP). Seizure onset was defined as the earliest occurrence of rhythmic or repetitive spikes in the cortex that was distinctive from the background activity, and that evolved in frequency and morphology (figure 2). The EEG onset of the seizure in the cortex was marked as 'unequivocal EEG onset' (UEO). SEEG segments were clipped to include 10 min before the UEO and 10 min post-termination of the seizure. Seizures were classified into focal aware seizures (FAS), focal impaired awareness seizures (FIAS), and focal to bilateral tonic-clonic seizures (FBTCS) [35]. A seizure was classified as electrographic seizures (ES) when it lacked any clinical correlate, and the ictal EEG pattern lasted longer than 10 s. A minimum of three seizures selected randomly among each seizure type (ES, FAS, FIAS, FBTCS) per patient, and overall 3–10 seizures per patient were included for analyses in the study. Multiple epochs of 5–10 min duration of interictal SEEG were selected and labeled as 'interictal state.' State of vigilance (SOV: sleep or awake) for the interictal epochs were determined based on video and visual SEEG inspection. No effort was made to classify sleep stages as the study was aimed to distinguish the thalamic ictal from interictal states. Similarly, no effort was made to determine the presence (or absence) of epileptiform spikes in the interictal epoch as the study was focused on determining the interictal state independent of presence (or absence) of a spike.

Zoom In Zoom Out Reset image size

Line length (LL) was used to detect the ictal activity objectively [36] (figure 1(B)). The LL feature was derived as a simplification of the running fractal dimension of a signal. It measures the length of the signal in a particular window and compares it to a threshold. The length of the signal is dependent on the amplitude and frequency of the signal, making this feature highly suitable to sense changes in amplitude and/or frequency that typically occur during seizures. The change in LL during the peri-ictal period was considered significant when the measured LL in that segment was greater than 2SD compared to the mean of a 3 min baseline segment for at least 10 s. Preprocessing steps involved detrending within 2 s, removing DC-drift (detrend), reconfiguring into a bipolar montage by linking adjacent contacts from the proximal to distal part of the electrode, and removing 60 Hz line noise with a notch filter (2nd order Butterworth zero phase shift filter—designfilt and filtfilt). Subsequently, the LL was calculated on this preprocessed data over 0.25 s windows with 50% overlap as shifting windows, and the resultant LL data was smoothed using the Matlab function 'movmean' over 2 s windows. The LL was calculated for each segment with m samples using the following formula:

$$\begin{equation}{\text{LL}} = \frac{1}{K}\mathop \sum \limits_{j = n - N}^n |{x_{j + 1}} - {x_j}|\end{equation} \tag{ 1 }$$

Where n is the number of sample points within 0.25 s, and the x_j is the sample point.

When LL detections were discontinuous,multiple close lying LL detections (<5 s between detections) were considered as a single seizure event. The algorithm detected 'seizure length' was the duration between LL detected onset (increased distance between preonset to onset) and LL detected offset (increased distance between termination and post-termination) (supplementary methods (available online at https://stacks.iop.org/JNE/17/066004/mmedia)).

The electrophysiological features of ANT involvement during a seizure are not well described. Thus, we have adopted a manual offline seizure detection approach based on LL features described above. LL was performed on all the cortical channels demonstrating visually determined UEO and early seizure spread in the cortex. Seizures in ANT were reported as 'not detected (ND)' when the specified changes in LL as described above, were absent in the first 120s after LL detected seizure in the cortex (LL UEO).

Two methods were used to measure the ictal changes of the thalamogram, (1) calculate the power-spectra with discrete wavelet technique with relative wavelet energy (DWT-RWE) and (2) multiscale spectral entropy (MSE) [26, 37, 38]. The time-frequency decomposition of field potentials was performed with the DWT—'db4', RWE, 6 levels. The DWT provides a non-redundant, highly efficient wavelet representation and direct estimation of local energies at the different relevant scales. The motivation for selecting DWT is based on prior studies suggesting this method can be an optimal tool for online seizure detection that can be translated into the implantable neural prosthesis [11, 39, 40]. The DWT-RWE was calculated on 4 s windows, with 3 s overlap, for 1–2 Hz, 2–4 Hz, 4–8 Hz, 8–16 Hz, 16–32 Hz, and 32–64 Hz from awake, sleep and seizure segments of the thalamic signal.

We calculated MSE to characterize the temporal predictability of a time series across several time scales in the thalamus, serving as an index of its capacity for processing information. Estimation of MSE was done in a windowed manner on short epochs of 4 s each, with 75% overlap after downsampling data to 200 Hz according to the original implementation of the algorithm and averaging across the 10 scales to indicate the mean MSE (https://archive.physionet.org/physiotools/mse/tutorials) [41].

We estimated MSE for one-dimensional discrete time series (x) by constructing consecutive coarse-grained time series, determined by the scale factor, τ, according to equations (2) and (3),

$$\begin{equation}y_j^{\left( \tau \right)} = \frac{1}{\tau }\mathop \sum \limits_{i = \left( {j - 1} \right)\tau + 1}^{j\tau } {x_i}\end{equation} \tag{ 2 }$$

where τ is the scaling factor, 1 ⩽ j ⩽ N/τ

$$\begin{equation}{\text{MSE}}\left( {m,r} \right) = - {\text{ln}}\frac{{{C^{m + 1}}\left( r \right)}}{{{C^m}\left( r \right)}}\end{equation} \tag{ 3 }$$

where ${C^m}\left( r \right) = \frac{1}{{\left( {N - m} \right)\,}}\mathop \sum \limits_i^{N - m} C_i^m\left( r \right)$, r is a coefficient of tolerance, m is the pattern length and N is the length of the data. Costa et al have shown that MSE had better statistical validity for m = 1 or 2 and r being between 0.1 and 0.25 [41]. In the present study we used m = 2, r = 0.15, 1 < τ ⩽ 10.

MSE was calculated across the length of the ANT SEEG signal on moving window epochs (awake, sleep, and seizure segments) like DWT analysis.

Seizures can induce changes in behavioral states (like arousal or altered vigilance) that can increase power at narrow frequency bands. To evaluate if the temporal and spectral features (MSE and DWT-RWE) can distinguish the thalamic ictal state from other interictal states of vigilance (awake and sleep), we used two machine learning algorithms—random forest (RF) and random kitchen sink (RKS) [42, 43]. The two classifiers RF and RKS are both suited for testing linear and nonlinear associations between the predictors and the response variables. One of the main reasons for choosing these algorithms is the nonlinear nature of brain dynamics [44, 45]. Evidence also states that linear methods are more suited for limited data and limited knowledge about the data available, while nonlinear methods are suitable to find a potentially more complex structure in large amounts of data. Specifically, when the source of the data to be classified is not well understood, it is preferred to use methods that are good at finding nonlinear transformations of the data [46]. In the current scenario, the recordings from the anterior nucleus of the thalamus and their electrophysiological behavior in human beings are being tested for the first time in SEEG. Also, in the future, when these methods are to be considered for building robust brain–computer interfaces for neuromodulation, the large data sets would be required to be trained using efficient models for implementation. The current nonlinear random RF models have a reliably higher accuracy among other algorithms, run efficiently on large databases, give estimates of what variables are important in the classification, has methods for balancing error in class population unbalanced data sets witnessed in real-time data evaluation, and can finally be extended to unlabeled data, leading to unsupervised clustering on new data sets (https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm). We further used the nonlinear RKS model as a validation tool, which by selecting many random non-linearities, can achieve the same accuracy in significantly less training time, a feature that can be exploited in real-time detections [43].

Seven features, namely: DWT-RWE in six frequency bands and the average of the first 10 scales of the MSE of awake, sleep, and 20 s of seizure segments (each seizure provided 16 data points) from the ANT were used as predictive features for RF and RKS machine learning algorithms. The seizure class contained 373 data points, and the interictal class had 2850 data points (approximately 4 h of data), each being represented by 7 features. The data was unbalanced, and hence the smaller data was split in 80%:20% for training and testing data. 300 data points were chosen at random for training, and 73 data points were chosen for testing purposes from each class. The RF classification technique was applied using all the features to calculate the accuracy, precision (positive predictive value), and recall (sensitivity) values. To ensure that training and testing were done in an unbiased manner, the process was repeated 500 times by choosing the training and test samples randomly during each iteration, and the average values of the accuracy, precision, and recall were used. The entire procedure used 5, 10, 15, and 20 decision trees, and four sets of results were presented that indicated that we have an improved classification with an increasing number of decision trees [47]. During the construction of the trees, Gini impurity was used to make decision about which variable to split at each node. For each feature, the sum of the Gini decrease across every tree in the forest is accumulated every time that feature is chosen to split a node. The sum is divided by the number of trees in the forest to give an average called the permutation feature importance.

As a validation of the results of the RF prediction model, we tested the same features using RKS learning. The RKS method involved a nonlinear mapping of the features on to a higher dimensional space and enabled the features to become linearly separable in the higher dimension. The feature mapping was done using a real Gaussian function as the radial basis function (RBF) kernel. This, in combination with a regularized least square algorithm for regression, allowed us to obtain a simple classifier that can be used for real-time applications. Similar to RF analysis, RKS was used to classify the given samples into two classes: seizure and interictal states (awake and sleep). RKS classification technique was applied using all the features to calculate the accuracy, precision (positive predictive value), and recall (sensitivity) values. The seven features (DWT-RWE in 6 bands and the average of the first 10 scales of MSE) of each data point were transformed into a 100 dimension feature space using random feature mapping. The mapping function performed both cosine and sine transformations on each data point, thus effectively converting each data point from a seven dimension vector into a 200-dimensional vector. Using the data in the new feature space, regression analysis, and regularized least square (RLS) loss function based classification was performed using an RBF based Gaussian kernel classifier with a hold-out cross-validation. The random mapping and the classification analysis was done using grand unified regularized least squares (GURLS) [48]. Similar to RF, the RKS was repeated 500 times by choosing the training and test samples randomly during each iteration, and the average values of the accuracy, precision, and recall were reported.

Since DWT-RWE showed non-normal distribution, we used the Wilcoxon rank-sum (WRS) method to test how the ictal spectral patterns differed (i.e. increase or decrease in DWT-RWE) from interictal (awake and sleep) state [49]. Five DWT-RWE samples of moving windows of 4 s length, with 4 s overlap, were input into the WRS test to form a time-related spectral power distribution subset from the seizures with ANT recruitment separated by seizure types for every frequency band. These subset distributions were compared to awake and sleep interictal distribution with the Wilcoxon test, and false discovery rate (FDR) correction was applied to the resulting p values [50].

Out of the 13 patients that had thalamic implants, 10 had accurate localization (77%) in the ANT and were included for subsequent analysis (table 1). None of the patients had any ischemia or hemorrhage in the subcortical regions, as evident from the CT brain performed after the explanation of the depth electrodes. The details of the accuracy and safety of research thalamic implantation at our center have been reported elsewhere [28].

Seventy-nine seizures from 10 subjects were analyzed in this study. Based on the clinical consensus among the epileptologists, the seizure onset zone was determined to be medial temporal (2), temporal pole (1), and temporal plus (2), with the plus representing additional seizure foci (orbitofrontal or insula) outside the medial temporal region. Five patients had multifocal epilepsies with seizure foci in bilateral medial temporal regions (2), medial temporal and parietal, or insular and superior temporal gyri (2). The distribution of seizure types was: ES (26), FAS (18), FIAS (29), and FBTCS (6).

The detection latencies of the seizures at the seizure onset zone channels were earlier (latency: 4.57 ± 10.28 s, range: −3.97–63.79 s) than in the ANT (latency: 16.89 ± 21.18 s, range: −0.99–92.98 s). The mean detection latencies in the ANT varied between seizure types: FAS (34 s), FIAS (12.04 s), and FBTCS (5 s).

Following seizure onset, the first 20 s in the ANT showed a progressive increase in spectral power in lower frequencies (4–16 Hz) for all seizure types (figure 1(C)). In addition, FIAS and FBTCS showed a sustained increase in 32–64 Hz. The peak changes in theta band (4–8 Hz) were seen around 26–31 s after seizure onset in the cortex. Also, we observed a decrease in MSE, implying that the degree of randomness of the thalamic neural activity reduced during the time of the seizures. The changes were noted for all types of seizures (ES, FAS, FIAS, and FBTCS).

Both RF and RKS were able to parse the thalamic seizure state from interictal states. For RF, the precision (positive predictive value) and recall (sensitivity) was 93.9% and 93.8%, respectively. For RKS, the precision and recall were 97.3% and 97.6%, respectively. In figure 3, the average values of the confusion matrices, and performance metrics that were performed over 500 iterations for both RF and RKS were highlighted. The performance was reliable across the 500 iterations for both RF and RKS with a narrow variation in the range of performance (accuracy of RF: 87%–99% and RKS: 84%–100%) (other performance measures are summarized in figures 4(A) and (B)). Permutation feature importance was done for the RF model in order to verify if the seven features (DWT-RWE and MSE) were relevant for the classification (figure 3(A)). We observed that all the features were relevant, with the theta and beta frequency bandwidths having the highest permutation importance. For the RKS model, the features were transformed by the nonlinear RBF kernel function, and the internally transformed features were used for classification. Since a nonlinear transformation was being performed, the feature relevance analysis was non-trivial in the transformed domain and was beyond the scope of this paper. To evaluate the significant differences between the features, we used univariate ANOVA that showed that all predictors differed significantly from each other except MSE and L2-4 spectral features (figure 3(B)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After establishing the relevance of the seven features in classifying seizures and the interictal data, the seven-dimensional data was visualized as a two-dimensional map using t-distributed stochastic neighbor embedding (t-SNE) algorithm. The 2-D map of the features was plotted using t-SNE with Mahalanobis distance metric (figure 3(D)). It can be observed that the DWT-RWE and MSE clusters during seizures were distinct from the interictal DWT-RWE and MSE, thus strengthening the probability of correct classification of thalamic states as seizure or inter-ictal. The ictal clusters were further noted to be sub-grouped based on the seizure types, i.e. ES, FAS, FIAS, and FBTCS, segregated into the distinct color clusters in the t-SNE plot (figure 3(D)). The receiver operating characteristic (ROC) curves show that both RF and RKS can discriminate the interictal and the ictal states very well (AUC of RF: 0.94 and RKS: 0.96) (figures 3(E) and (F)).

Overall, these results demonstrate that features extracted in the time and spectral domain from the thalamic LFPs can be used to classify seizures and interictal states in focal epilepsy.

Although multiple spectral changes characterized the thalamic ictal state, the emergence of theta band was significantly prominent and sustained for seizures that remained focal. Our findings are in agreement with recent studies that demonstrated increased power in lower frequencies in the ANT at the seizure onset [57]. Vertes et al demonstrated the presence of rhythmic burst firing of neurons in the ANT that was synchronous with hippocampal theta rhythm, and the rhythm resonated within the Papez circuit (including medial septum hippocampus and ANT) [58]. Prior studies have demonstrated that theta rhythm is associated with increased seizure threshold, is anti ictogenic, and synchronization in theta rhythm can ameliorate epileptic discharges [59, 60]. Further work is necessary to establish the role of theta rhythm in seizure propagation.

With the increase in multifocal and difficult-to-localize epilepsies, there is a need to develop effective neuromodulation therapies as surgical resection is often unavailable. The ANT mediates cortical-subcortical interactions between the limbic system and the brainstem. These pathways are associated with a generalization or bihemispheric propagation of convulsive seizures [61]. Thus the thalamic hub is an attractive 'choke point' for modulating the distributed epileptogenic network and disrupting the seizure genesis [62, 63]. Indeed stimulation of the ANT confirmed modulation of the limbic system. Furthermore, recent case studies highlighted the feasibility of the RNS system in modulating the thalamocortical nodes during seizures [53, 64]. Our work is a step forward in identifying reliable seizure markers in the thalamus. The ability to predict seizure types at the onset might allow developing stimulation parameters that are specific for seizure types.

In the present study, the data was trained on seizure clippings as discrete events, and the event detection was considered as a classifier and not as a prediction. Future studies should validate the findings on continuous extended periods of EEG recordings from the thalamic subnuclei to affirm that such algorithms can be implemented in real-time usage. Secondly, the features used for machine learning were limited to entropy and spectral measures. Other features like the synchronization patterns between the seizure onset zone and the thalamus are likely to influence the detection rate of seizures. Lastly, we noted that these machine learning tools were time expensive, where our data sets were trained and tested over hours. For neuromodulation based on online event detection, the algorithms need to perform faster while retaining the higher accuracy and sensitivity.