For the current clinical study, we analyzed 42 subjects: 14 patients (9 male) suffering from frontal focal epilepsy, 14 patients (5 male) suffering from idiopathic generalized epilepsy (they all meet the criteria Juvenile Myoclonic Epilepsy or presumed genetic, as the new terminology stands; ILAE, 1989), and 14 (7 males) healthy subjects. The demographics of the subjects is summarized in Table 1.

All patients were classified independently by two neurologists. For the diagnosis and classification of seizures and epileptic syndromes, the clinical and electroencephalographic classifications of the ILAE were applied (ILAE, 1989). Patients were seen consecutively during the period between May 2009 and December 2011 in the outpatient clinic of neurology at the University General Hospital of Ciudad Real, with an average evolution track of 10 years, and were free from epileptic crisis during the last 6 months prior to the recording. They were all free from mental retardation, connate anoxia history, history of head trauma or meningitis in early childhood. Each clinical history gathered epidemiological data (sex, age of onset of the first seizure, personal history, pregnancy, and delivery), clinical data (type of seizures and neurological findings), and treatments used. These data were combined with additional examinations [EEG, neuroimaging studies, magnetic resonance imaging (MRI) of 1.5 or 3 Tesla following the protocol for epilepsy that includes 3D-T1 isotropic voxel sequences of 1 and 2 mm axial and coronal slices with T2, FLAIR, and T2^* sequences], to establish if the etiology was idiopathic, symptomatic, or cryptogenic. The entire cohort of patients was also evaluated using a comprehensive neuropsychological test battery. The patients took a range of medications with different mechanisms of action. In this study, patients were chosen in such a way that the medication was similar and controlled between both epilepsy groups (Niso et al., 2015). All data were analyzed anonymously and ethical approval was granted by the Local Ethics Committee of the Teaching General Hospital of Ciudad Real.

For all subjects, we recorded brain activity during resting state with MEG at the Center for Biomedical Technology, CTB, Madrid. MEG recordings were obtained using a 306-channel whole head Elekta Neuromag MEG system (Elekta Oy, Helsinki, Finland), comprising 102 magnetometers and 204 planar gradiometers in a helmet-shaped array covering the entire scalp, while subjects were seated inside a magnetically shielded room (Vacuumschmelze GmbH, Hanau, Germany). Eye movements were monitored by simultaneously recording the electrooculogram with three Ag/Cl electrodes, two above and below the right eye and one at the right earlobe (ground reference). Four head position indicator (HPI) coils, whose data were used to correct head movement during the session, were placed on the scalp, appropriately spaced in the region covered by the MEG helmet. The locations of the nasion, two pre-auricular points, and the 4 HPI coils were digitized prior to each MEG study using a 3D-digitizer (FASTRACK; Polhemus, Colchester, VT) to define the subject-specific Cartesian head coordinate system. One hundred to two hundred additional anatomical points were digitized on the head surface to provide a more accurate shape of the subject's head.

MEG data were acquired (sampling rate of 1 kHz, on-line band-pass filter ranging from 0.10 to 330 Hz) during a single session of resting state (Niso et al., 2015). For each session, we analyzed data obtained during 10 min with eyes closed. These data were free from blinking artifacts and interictal epileptiform discharges. External noise was removed, as recently suggested in MEG literature (Gross et al., 2013), by using the temporal extension of Signal-Space Separation (tSSS) (Taulu et al., 2005) implemented with the MaxFilter software (version 2.0 ElektaNeuromag; sliding window of 10 s, subspace correlation limit of 0.98; Hillebrand et al., 2013), and configured with an inner expansion order of 8. We also applied movement compensation to our data, to correct for subject's head movement during the recording. Spikes and spike-wave complexes were observed in the data but were removed from the analysis since our aim is to classify subjects during the absence of epileptiform activity. This is a key aspect of our study since it is not always easy to record epileptiform activity during short periods of recording times. Following this procedure, and within the 10 min with eyes closed, we selected the 40 most stationary segments of 5 s length for each patient by using the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test for stationarity (Kipiński et al., 2011).

For each temporal segment, we consider the information of the 102 magnetometers in the MEG system, i.e., each segment of 5 s contains information of 102 different channels.

We compute four different features for each segment, namely the total and relative power spectral densities (PSD), the phase-locking value and the phase-lag index, described below. These features are calculated for every segment and sensor (magnetometer) of each subject. All these measures are calculated for a range of frequencies of interest in the interval (4–40) Hz, in steps of 2 Hz. As a result, the computation of each feature leads to a matrix of dimensions 102 x 18 (sensors x frequency intervals), which is used as input to the machine learning algorithm discussed in Section 2.4 for further classification of the MEG segments and subjects.

We compute total and relative PSD of the MEG data (40 segments for each patient) for the range of frequencies of interest. Relative PSD refers to the total PSD on each sensor on a certain frequency band divided by the total PSD on that sensor. In Figure 2, we show an example of the relative PSD for one MEG segment of a given subject. This image contains the spectra of the 102 sensors (y-axis) in 18 frequency intervals of 2 Hz (x-axis).

We analyze synchronization and functional brain connectivity from the MEG data by computing two bivariate phase synchronization measures: the phase-locking value (PLV) (Lachaux et al., 1999; Mormann, 2000) and the phase-lag index (PLI) (Stam et al., 2007). The use of PLV and PLI as features is motivated by the fact that epilepsy has been historically seen as a functional brain disorder associated with excessive synchronization of large neuronal populations (Jiruska et al., 2013). Increasing evidence shows that epileptiform phenomena, particularly seizures, follow a characteristic dynamical evolution of synchronization. Desynchronization is often observed preceding seizures while high levels of synchronization are observed toward the end of seizures. For computations, we use the HERMES toolbox (Niso et al., 2013) (freely available at http://hermes.ctb.upm.es/). Amplitude-based measures such as mutual information or measures of generalized synchronization are less suitable to analyze synchronization regimes in which the amplitudes remain only weakly correlated. Instead, PLV is able to detect weak synchronization regimes, where the phases of the oscillatory component are coupled but the amplitude may not be (Rosenblum et al., 1996; Hramov et al., 2005). In turn, PLI is a measure of asymmetry for the distribution of phase differences between two signals (Stam et al., 2007). PLV is sensitive to both zero and non-zero lag synchronization, while PLI is only sensitive to the latter one. Zero lag synchronization is usually regarded as the result of the same neural source being picked up by two different sensors. In addition, zero lag synchronization can be achieved if two neural sources are indirectly connected through a third one, acting as a dynamical relay (Vicente et al., 2008), which would be overlooked should we only focus on PLI. Therefore, we include both PLV and PLI as features in our study so that we can elucidate whether information from only direct connections (PLI) is enough for a good discrimination or, instead, including information on both direct and zero lag functional connectivity is best for this purpose.

We compute the PLV and PLI on bands centered in the frequencies ranging from 4 to 40 Hz with a resolution of 2 Hz, which represents a reasonable trade-off between two mutually opposed requirements. First, the frequency intervals should be narrow enough to allow for proper phase reconstruction (Thiel et al., 2006). Second, selecting extremely narrow intervals would require a very high-order filter, which may distort part of the data segments due to border effects. It is common practice to apply a threshold of significance to the values of PLV and PLI (Kim et al., 2013). Here, we keep the 40% (50%) most significant values of PLV (PLI) and zero the rest, keeping a fixed proportion of the highest links. However, we have checked that the results do not change qualitatively by choosing other thresholding levels. We use the average PLV and PLI as estimators of the functional connectivity at each sensor location, within the frequency bands mentioned above. PLV and PLI are computed for each pair of sensors, such that the PLV_{i, j} (PLI_{i, j}) is the PLV (PLI) between sensors i and j. The corresponding average PLV (PLI) of sensor i is then the average of the PLV_{i, j} (PLI_{i, j}) sum over all j, with i ≠ j. In Figure 3, we show the computed average PLV for a single MEG segment of a given subject. Similarly to the procedure described for the total and relative PSD, the results of the PLV or the PLI are passed in matrix format of dimensions 102 x 18 (sensors x frequency intervals) as input to the machine learning algorithm for the classification of the MEG segments and subjects.

Here, we use randomized neural networks as the classification algorithm. In particular, we employ ELMs, which are supervised machine learning models based on a random mapping of the inputs into the high dimensional space of the neural network. ELMs aim at learning the conditional expected value of a m-dimensional output y(n) as a function of a d-dimensional input x(n), based on a number of examples n = 1, …, N. The particularity of this approach is the use of an intermediate hidden layer with D neurons, where inputs are randomly mapped into a nonlinear form, yielding a new, transformed D-dimensional r(n) space. The elements of the hidden layer are often denoted as neurons since these approaches were initially developed in the framework of neural networks.

The ELM model can be written in the following general form (Huang et al., 2006; Ortín et al., 2015): (1)r(n)=F(ηWinx(n)), (2)o(n)=Woutr(n), where r(n) and o(n) are the states of the hidden layer neurons and the predicted output values, respectively. W_in is the d × D random input weight matrix that maps the input (dimension d) onto the hidden layer (dimension D), F is the hidden layer nonlinear activation function, and η is an input scaling factor. We use F(z) = sin²(z + ϕ), where z is the mapped input signal, as the nonlinear activation function. For this activation function, the operating point and hence the degree of nonlinearity can be adapted by varying the phase ϕ. Machine learning models with this nonlinear function can be easily implemented in optoelectronic hardware, allowing for ultra-fast processing speeds (Larger et al., 2012; Soriano et al., 2015a).

Learning from data is efficiently achieved through the random mapping technique, since the only weights trained in this approach are those corresponding to the output connections, W_out. W_out are usually computed by minimizing the squared error between the real y(n) and the predicted o(n) output values. In our case, we have used a simple linear regression (Huang et al., 2006). The input data x(n) are first multiplied by W_in and then normalized to zero mean and unity standard deviation. In our study, the input data correspond to the flattened matrices of the computed features. Thus, the input features for each segment typically have a dimension d = 1, 836 (102 x 18).

The values of W_in are continuous, drawn from a uniform distribution over the interval [1, −1]. Since they are randomly chosen, the performance of a single realization of the algorithm might depend on the particular values of these random input weights. To avoid this uncertainty, we use an ensemble classifier that produces statistically robust results. This ensemble classifier is a combination of the results from a pool of independent realizations of W_in. Here, we use ensemble classifiers that combine the results of 100 different realizations using a simple majority rule.

One of the main advantages of this classification algorithm based on random mappings is that the parameters of the algorithm can be easily retrained if new data is incorporated, either from new subjects or more measurements from the same subjects. This property is particularly suited for applications where data is acquired sequentially, e.g., new subjects being incorporated to the study.

The performance of the classifiers is measured in terms of four quantities defined as the number of true positives (TP), the number of true negatives (TN), the number of false negatives (FN), and the number of false positives (FP). These quantities are often displayed as a confusion matrix, following a table layout. As shown in Table 2, the columns of the confusion matrix represent the instances in a predicted class while the rows represent the instances in a true class.

From the confusion matrix, it is easy to compute the total proportion of correct predictions (accuracy), the percentage of positives that are correctly identified as having the condition (sensitivity), and the percentage of healthy subjects that are correctly identified as not having the condition (specificity), defined as follows: (3)Accuracy = TP+TNTP+TN+FP+FN , (4)Sensitivity = TPTP+FN , (5)Specificity = TNFP+TN . At the subject level, the results are presented in the form of the confusion matrix. At the segment level, we make use of the receiver operating characteristic (ROC) curves to present the results. In statistics, the ROC curve is used as a graphical tool to illustrate the performance of a binary classifier as its discrimination threshold is varied. The ROC curve is created by plotting the probability of detection (sensitivity) against the probability of false alarm (1-specificity) at various threshold settings. Thus, each point on the ROC curve represents a sensitivity/specificity pair corresponding to a given discrimination threshold. Performance is measured by the area under the ROC curve. This area under the curve (AUC) measures the ability of the classifier to correctly identify those subjects with and without the epileptic condition. An area of 1 represents a perfect classifier, while an area of 0.5 represents the same classification as in a random guess.

For the training and test of the ELM classifiers, we follow a leave-one-out cross-validation procedure. This procedure is a standard way to evaluate the performance of machine learning algorithms, especially when the number of subjects is small as in our case. Leave-one-out cross-validation involves using one subject as the test set for the classifier trained using the remaining subjects. In our case, we have 14 subjects for each condition (healthy, focal, generalized) that we divide in 14-folds. In this manner, one subject of each condition is tested in every fold and all subjects are tested exactly once after the leave-one-out cross-validation. In this way, we do not overestimate the accuracy of the algorithm using intra-patient information to classify the segments. The subject-wise approach, in contrast to a segment-wise one, estimates better the accuracy of the algorithm (Saeb et al., 2016).

The ELM classification algorithm is first trained to classify MEG segments using as input one or more of the features defined in Section 2.3. After the training and for each fold, the algorithm assigns a condition to each segment of the test subject. It then decides on the condition of the test subject by using a majority rule over the results for the segments of this test subject. This process is repeated for each of the 14 test folds and the final accuracy of the classifier is obtained by averaging over the individual fold accuracies. We note that the training and test sets for each fold contain records from non-overlapping subjects since our goal is to build a classifier that can generalize to identify epilepsy traits rather than the subject itself.

The parameters of the ELM algorithm are optimized for the identification of epileptic and healthy segments. Since the current study is performed on a total of 42 subjects, with 40 temporal segments each, this amounts to 1,680 (42 x 40) examples. This reduced amount of training examples limits the number of neurons that can be used in the hidden layer of the ELM, which must be much smaller than the number of examples. In this work, we present results for a hidden layer with D = 350 neurons since we have checked that it gives a consistent performance.

The other parameters of the ELM algorithm described in Section 2.4 that need to be tuned are the input scaling η and the phase ϕ. To do so, we test the accuracy of a binary classifier to identify the condition of epilepsy in the MEG segments as a function of the algorithm parameters. Figure 4 illustrates the results for a η and ϕ parameter scan. From the figure, it is clear that the parameter values η = 0.1 and ϕ = 2.1 offer a good classification accuracy and will be kept fixed throughout the rest of the manuscript.

Here, we present the results for the first stage of the automated classifier, trained to distinguish between MEG segments of healthy and epileptic subjects (without specifying the type of epilepsy). We train an ELM independently for each feature to determine the best discriminant feature between healthy and epileptic subjects, each time taking only one of the calculated features as input. Figure 6 shows the ROC curves obtained for the identification of the MEG segments of epileptic subjects for each feature. We find that the best results are obtained when we use the relative PSD as input to the ELM classifier. In turn, the total PSD also performs better than the phase synchronization measures. Finally, PLV is superior to PLI, which yields classification values close to chance.

The ROC curves of Figure 6 illustrate the algorithm performance for the classification of segments. It is, however, also interesting to discuss the results in terms of subjects. To this end, we check the condition (healthy or epileptic) of the corresponding segments of each subject and assign a condition to the subject following a majority rule vote over his/her segments (see Section 2.6 for details). We do it only for the relative PSD, as it gives the best results at the level of segments. The test performance of the classifier at the subject level is summarized in the confusion matrix shown in Table 3. This confusion matrix has an accuracy of 0.9, a specificity of 0.86, and a sensitivity of 0.93. These results show that it is possible to identify epileptic subjects from interictal MEG resting activity with high accuracy using an ELM classifier. In the test phase of the cross-validation procedure, the ELM classification algorithm is able to properly classify 38 out of the 42 subjects. The nonlinear transformation carried out by the ELM is essential to obtain these results, as performing linear regression over the original input features, without applying the ELM algorithm (results not shown), yields classification accuracies close to chance.

We now evaluate which features obtained from the MEG time-series are able to best distinguish between generalized and focal traits. To this end, we train a classifier with only the epileptic subjects. We find that the relative PSD again outperforms the other features for the identification of the epilepsy type, while the PLV is the synchronization measure with the best performance. The overall classification accuracy is, however, low even for the relative PSD. In order to try improving the classification of epilepsy types, we evaluate the performance of the ELM classifier when we restrict the frequency content of the features to specific frequency bands (Niso et al., 2015). In particular, we focus on the theta (4–8 Hz), alpha (8–14 Hz), beta 1 (14–20 Hz), beta 2 (20–28 Hz), and low gamma (28–40 Hz) bands. We find that the relative PSD, using only the beta 1 frequency band, helps best to distinguish between generalized and focal epilepsy. This means that the ELM algorithm classifies better the types of epilepsy if the features are restricted to frequency bands before automated classification. We then search if there are pair-wise combinations of (frequency restricted) features that can further improve the classification. The best classification performance, as measured by the AUC (Table 4), is obtained by using a suitable pair-wise combination of the relative PSD (beta 1) and PLV (beta 2) of the frequency restricted features. The choice of pair-wise combinations is motivated by the hypothesis that these two different features contain complementary information of the underlying data. Indeed, PSD reflects the local activation of each sensor data, whereas the average PLV reflects the statistical relationship between the MEG activity in this sensor and that in the rest of the sensors. Namely, the former is an estimation of local activity, whereas the latter estimates changes in global connectivity.

For the combination yielding the highest AUC in Table 4 (beta 1 band of the relative PSD and beta 2 band of the PLV), the corresponding ROC curve for the MEG segments is presented in Figure 7. In addition to the ROC curve for the optimum combination of input features, we present the ROC curves for these input features considered individually. By comparing these ROC curves, it is apparent that the AUC of this combination is higher than the AUC for the beta 1 band of the relative PSD, alone.

Similarly to the procedure described in Section 3.1, the test results at the subject level for the classifier trained to identify the type of epilepsy are presented in Table 5. These results correspond to the optimal combination detected above of the beta 1 relative PSD and beta 2 PLV. The accuracy of this binary classifier is 0.93 (specificity of 0.86 and sensitivity of 1). We find that, in the test phase, the classifier is able to properly identify the type of epilepsy for 26 out of the 28 epileptic subjects studied.

As reported in previous sections, it is possible to distinguish between epileptic and healthy subjects (stage 1) and epileptic subjects with general or focal epilepsy (stage 2) with high accuracy. In this section, we report the test results for the two-stage classification procedure illustrated in Figure 5. This procedure allows identification of the three different conditions (healthy, generalized, and focal). The complete classification procedure contains two binary ELM classifiers. Thus, each test subject is first classified as healthy or epileptic. In the latter, a second classifier further discriminates the type of epilepsy.

We present the confusion matrix for the full classifier in Table 6. As input features, we choose the relative PSD for the first stage and the combination of beta 1 relative PSD and beta 2 PLV for the second stage. In the test phase of the final classifier, 36 out of the 42 subjects are properly classified. The accuracy of the full classifier for the three conditions is 0.86, with a sensitivity of 0.88, and a specificity of 0.86.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.