Automated detection of epileptic ripples in MEG using beamformer-based virtual sensors

Nine patients diagnosed with intractable epilepsy were selected for this study. All of them had simultaneous MEG and iEEG recordings. The subjects selected for this study presented interictal activity coming from one focal generator. Table 1 shows a summary of the available clinical information, seizure type, the area where the ictal electrode was identified by clinicians and the iEEG contact number.

MEG signals were acquired in a magnetically shielded room over 10 min with eyes closed, with the patient in a supine position in order to remove head movements. A whole-head 148-channel magnetometer system (4D-Neuroimaging/BTi San Diego, CA) sampling at 678.19 Hz (bandwidth DC to 250 Hz) was used to acquire MEG signals and, simultaneously, iEEG was recorded using the same sampling rate and bandwidth conditions. For each patient, one depth electrode was implanted in the target region. Each electrode was 1.1 mm in diameter and included six or eight leads, 2.4 mm in length, 5 mm apart (AD-TECH, Racine, WI). The relative position of the patient's head in respect to the MEG sensors was recorded continuously using head-localization coils. Signals were imported into MATLAB using the Fieldtrip toolbox [29]. MEG signals were filtered to remove ocular, cardiac and metallic artifacts using automatic blind source separation algorithms [30, 31], and visually inspected to select 200 s of interictal periods for posterior analysis. These algorithms have proven to be effective in the localization of epileptic sources [32–34]. The general recommendation when visualizing HFOs is to use a sampling frequency that at least five times the frequency of the oscillation of interest [8] and, for this reason, the maximum analysis frequency was set to 120 Hz. MEG and iEEG signals were band-pass filtered from 1 to 120 Hz to obtain the whole spectrum (MEG_WS) within these two frequencies and from 80 to 120 Hz to obtain the high frequency spectrum (MEG_HF) in both cases using a finite impulse response filter (80th order).

The patient's index points and head shape were digitized with 3Space Fasttrack (Polhemus, Vermont, USA) prior to each measurement. The nasion, an anatomical landmark, and the left and right ear canal points served as index points and were used to define a right-handed coordinate system, called the headframe coordinate system: the x-axis points to the front, the y-axis to the left and the z-axis to the top of the head.

The scalp and brain meshes for each subject were obtained by aligning and warping the default anatomy provided by the Montreal Neurological Institute with Brainstorm software [35]. The volume conduction model was obtained with the single-shell algorithm provided by Fieldtrip in which a realistic single-shell model based on lead field expansion is computed [36].

The proposed algorithm is divided into two main steps:

• In the first step, beamforming was applied to the whole head within a coarse grid of MEG-VS to delimitate the ROI. As beamformer computation time increases with higher resolution, this step permitted delimitating the area of analysis efficiently.
• In the second step, a second beamformer in the ROI was computed and a finer MEG-VS grid was obtained. The automatic detection of ripples was performed with a modified version of the algorithm developed by [10]. This algorithm uses the Stockwell transform [29], which provides excellent time-frequency decomposition and reduces computational run time.

In order to assess the performance of the proposed procedure, ripples detected in simultaneous iEEG signals were used as the gold standard. To facilitate understanding of the whole procedure, a scheme is shown in figure 1 and a detailed explanation is provided in the next sections.

Zoom In Zoom Out Reset image size

2.3.1. First stage: whole-head beamformer and detection of the ROI.

2.3.1.1. Whole-head beamformer.

The first stage was performed to define the area where the epileptic high-frequency activity was taking place. To obtain this area, a linearly constrained minimum variance beamformer (LCMV-Beamformer) [37] was computed to obtain the weights of the spatial filtering using an evenly spaced (1.5 cm) grid of VSs including the whole brain volume of each subject. Beamformer weights (BF_w) were computed from MEG_HF using sliding windows of 100 ms with a 50% overlapping. For each time segment, virtual sensors (VS_HF) were obtained as the product of the MEG_HF signals with the filter weights projecting the time series along a unique dipole direction corresponding to the maximum variance, following equation (1):

$$\begin{eqnarray}&&\begin{array}{c}\text{V}{{\text{S}}_{\text{HF}}}\\ n\times s \end{array}=\begin{array}{c} \text{B}{{\text{F}}_{\text{w}}} \\ n\times m \end{array}\cdot \begin{array}{c}\text{ME}{{\text{G}}_{\text{HF}}}\\ m\times s\end{array}\end{eqnarray} \tag{ 1 }$$

where n is the number of VS_HF, m is the number of MEG channels and s is the length of the signal. The noise covariance was computed with the first 10 s of the unfiltered data [24]. Finally, MEG-VS signals were reconstructed by averaging the overlapped and successive time segments.

2.3.1.2. Detection of events of interest (EOIs).

For each computed VS_HF, EOIs were detected using an adapted version of the first phase of the algorithm developed in [10] for iEEG signals. For each VS_HF:

• (1)  
  The envelope of the band-pass signal was computed with the Hilbert transform.
• (2)  
  The mean and the standard deviation (SD) of the envelope was computed to set a threshold equal to mean  +  3  ×  SD.
• (3)  
  An event was detected when the envelope exceeded the threshold. The duration of the event was defined as the interval between the upward and downward crossings through half the detection threshold, that was 0.5  ×  (mean  +  3  ×  SD).
• (4)  
  If the duration of the event exceeded 6 ms, the event was qualified as an EOI.
• (5)  
  EOIs with an inter-event interval of less than 10 ms were merged.
• (6)  
  Events without a minimum of four peaks greater than 2  ×  SD from the mean baseline were rejected.

2.3.1.3. ROI selection.

By using the information provided by the EOIs, where a preliminary selection of high-frequency activity was performed discarding a fair amount of noisy events, a region with greater epileptogenic high-frequency activity could be delimitated. To do so, the number of EOIs appearing in each VS_HF was computed, all possible 3  ×  3  ×  3 subgrids were considered and the one showing the highest number of EOIs was selected as the ROI.

2.3.2. Second stage: detection of ripples inside the ROI.

2.3.2.1. Beamforming inside the ROI.

A new and finer grid (0.75 cm between VS_HF) was computed inside the selected cubic ROI, creating a volume of 5  ×  5  ×  5 points where new beamformer weights were calculated and time-series signals were obtained for each VS_HF by using the same procedure applied in the whole-head approach.

2.3.2.2. Detection of EOIs.

EOI detection was performed for VS_HF using the same algorithm explained in section 2.3.1. However, visual analysis of the selected segments revealed that the detection threshold of the algorithm was very restrictive and some of the events detected in iEEG that were observable in VS_HF were not being selected as EOIs due to a lower SNR and, therefore, failed to detect a considerable number of visible events in MEG-VS. Consequently, to perform the selection of EOIs and the subsequent detection of ripples, several thresholds' mean  +  factor_SD  ×  SD were evaluated by varying the factor_SD from 1.5 to 3 in steps of 0.25.

2.3.2.3. Recognition of ripples among EOIs.

The recognition of ripples among the selected EOIs was performed using a similar algorithm to the second step explained in [10], which discards EOIs elicited by activity other than oscillations [5, 13, 38]. The assumption taken for this step was to define a HFO as a short-lived event with an isolated spectral peak at a distinct frequency [39, 40]. The beamformer was computed from MEG_WS for each detected EOI  ±  0.25 s and VS_WS was obtained as the product of the MEG_WS with the filter weights projecting the time series along a unique dipole direction corresponding to the maximum variance as explained in section 2.3.1. The time-frequency response of the VS_WS was performed by using the Stockwell transform. An event was considered a ripple if its frequency response presented a clear distinguishable peak above 80 Hz. To measure this, the following steps were taken:

• (1)  
  The instantaneous power spectrum density (PSD) was obtained for the maximum point of the envelope (obtained by means of the Hilbert transform). An example of the computed time-frequency response and the corresponding PSD is given in figure 2.
• (2)  
  From the PSD, three parameters were computed (see figure 2):

  • (a)  
    The high-frequency peak (HiFP) was selected as the spectral peak of a possible HFO. The minimum frequency of this peak is 80 Hz and the maximum is 120 Hz.
  • (b)  
    The trough was defined as the minimum point in the range between 40 Hz and HiFP.
  • (c)  
    The low frequency peak (LoFP) was defined as the closest local maximum below the trough.
• (3)  
  The EOI under analysis was selected as a ripple if the following conditions were fulfilled:

  $$\begin{eqnarray}&&\frac{\text{Power}\left(\text{Trough}\right)}{\text{Power}\left(\text{HiFP}\right)}<0.8\end{eqnarray}$$

  $$\begin{eqnarray}&&\frac{\text{Power}\left(\text{HiFP}\right)}{\text{Power}\left(\text{LoFP}\right)}>0.5.\end{eqnarray}$$

Zoom In Zoom Out Reset image size

2.3.2.4. Whole grid analysis.

An additional factor was taken into account to discard ripples that appeared in just a few VSs and were isolated from their adjacent channels. All the detected ripples occurring at the same time segment (with a minimum of 50% of overlapping) were assigned to the same specific ripple event. The overlapping conditions are shown in figure 3. The relative power of the ripple in one virtual sensor RP(VS_i) was computed as follows:

$$\begin{eqnarray}&&\text{RP}\left(\text{V}{{\text{S}}_{i}}\right)=\frac{{{P}_{\text{event}}}\left(\text{V}{{\text{S}}_{i}}\right)-{{P}_{200\,\text{s}}}\left(\text{V}{{\text{S}}_{i}}\right)}{{{P}_{200\,\text{s}}}\left(\text{V}{{\text{S}}_{i}}\right)}\end{eqnarray}$$

where P_event is the power in the 80–120 Hz band of the ripple event in the ith VS and P_{200 s} is the power of the whole 200 s signal in the same frequency band in the ith VS.

Zoom In Zoom Out Reset image size

The energy of an event (E_r) was then computed as the sum of the maximum relative power of each VS belonging to the ripple with respect to the sum of the relative power of all the VSs in that same event:

$$\begin{eqnarray}&&{{E}_{\text{r}}}=100\cdot \frac{\underset{{}}{\overset{{}}{{\sum\nolimits_{n=1}^{n\text{VS}}}}}\text{RP}\left(\text{V}{{\text{S}}_{n}}\right)}{\underset{{}}{\overset{{}}{{\sum\nolimits_{k=1}^{125}}}}\text{RP}\left(\text{V}{{\text{S}}_{k}}\right)}\end{eqnarray}$$

where nVS is the number of VSs with a detected ripple. Once the E_r was calculated, the detected ripple was discarded if its value was not high enough. The energy threshold $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$ was defined to discard isolated events or those with low energy and, thereby, obtain better algorithm performance. Several values for $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$ were considered from 0% to 100% in steps of 1%.

Evaluation of the performance of the current algorithm was performed using iEEG signals as the gold standard. iEEG ripples were automatically detected using the time-frequency analysis algorithm described in [10]. For each pair of values considered in both thresholds (factor_SD and $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$) and for each automatically detected ripple in iEEG and MEG-VS, a co-occurring event was defined if these aforementioned ripples overlapped in time (see figure 3). Each ripple event in an iEEG was considered independently for each channel. In MEG-VS, all the VSs containing an overlapped event (figure 3) were considered as one individual event. This co-occurring ripple was considered a true positive (TP). An FP was defined as a detected ripple in MEG-VS that was not detected in iEEG. Finally, a false negative (FN) was defined as a detected ripple in iEEG but not visible in MEG-VS. The recall (R), precision (P) and F1 score (F1) were computed as:

$$\begin{eqnarray}&&R=100\cdot \frac{\text{TP}}{\text{TP}+\text{FN}}\end{eqnarray}$$

$$\begin{eqnarray}&&P=100\cdot \frac{\text{TP}}{\text{TP}+\text{FP}}\end{eqnarray}$$

$$\begin{eqnarray}&&\text{F}1=2\cdot \frac{P\cdot R}{P+R}.\end{eqnarray}$$

The F1 score measures the test accuracy and can be interpreted as a weighted average of the precision and the recall. The influence of both thresholds on the discrimination of ripples from noise was evaluated by observing increases of recall or precision. To select the optimum threshold values, a leave-one-out procedure was used. The values of energy and SD threshold were computed as a function nine times, each one leaving one of the subjects out. For each curve, E_r and SD values that produced higher F1 scores were selected. The values selected for the classifier, factor_SD and $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$ were selected using the median values for each of the leave-one-out functions.

As mentioned before, the ROI was selected as an area with a higher number of EOIs, and where the VSs belonging to the ROI were expected to be more active. In order to test this hypothesis, once the automatic selection was performed and the remaining noisy events were discarded, a comparison between the number of detected ripples inside and outside the ROI was performed. For each patient, the contralateral hemisphere was selected as the area outside the ROI and the same automatic procedure to detect ripples was applied. Results were compared between the two regions.

To evaluate the selected ROIs resulting from the first stage in the volume reconstruction model, the coordinates of the central virtual sensor were normalized into MNI space using Brainstorm software (figure 4) and the gyrus and lobes were obtained for this position in the normalized atlas (table 2). The central VSs of the nine subjects matched the lobe where seizures had been generated (table 1).

Zoom In Zoom Out Reset image size

To select the values of factor_SD and $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$ as the trade-off between precision and recall, the detection in iEEG and MEG-VS was compared when varying these two parameters, selecting the values that produced higher F1 scores, thereby maximizing both precision and recall. To obtain these values, a leave-one-out cross-validation approach was used, in which each subject was left out and the selected values were calculated with the rest. The results of the selected values for factor_SD and $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$ are shown in table 3. The obtained median values were used as the final selected values: 2.25 for factor_SD and 17.75% for $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$. The averaged precision/recall and F1 score curves for all subjects are shown in figure 5. A decrease of factor_SD permitted detection of a higher number of ripples, increasing recall but with the risk of introducing FP events (decreasing precision) with less energy that would appear simultaneously in a lower number of VSs. This is why a similar effect was found with $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$: higher precision but lower recall with higher values.

Zoom In Zoom Out Reset image size

The automatic detection was applied to the nine subjects using the selected threshold values. The values of recall and precision were 68.9  ±  0.8 and 79.2  ±  1.5, respectively (mean  ±  SD of all subjects). The resulting F1 score was 73.7  ±  0.6 (TP: 25.2  ±  10.0; FN: 11.4  ±  4.7; FP: 6.6  ±  2.5). The average number of detected events was 31.8  ±  12.5, which produced a rate of ripples per minute of 9.5  ±  3.7, similar to the value of 6.7  ±  4.8 obtained by visual inspection in [24] with MEG. The number of detected ripples in iEEG was 36.7  ±  14.7, producing a rate of 11.0  ±  4.4 ripples per minute. The detector was evaluated to ensure better-than-chance performance [41] using the Poisson random probability distribution proposed in [42]. Precision, recall and F1 score were significantly lower (p-value  <  0.0001) for the random distribution.

Figure 6 shows, as an example, MEG-VS time-domain signals for patient 1 after applying the beamforming algorithm. The same time segments of the iEEG channel with a higher ripple rate and the nine closest scalp-MEG channels are also shown. Ripples that are unrecognizable in scalp-MEG are observed in MEG-VS, and occur simultaneously in iEEG signals. Furthermore, the different possible situations (TP, FP, FN) of automatic detection of ripples in MEG and iEEG are shown. Table 2 summarizes the performance of the automatic algorithm for each patient along with the global performance using the aforementioned thresholds.

Zoom In Zoom Out Reset image size

Note that the automated detection used for iEEG signals was based in the automatic algorithm developed in [10], where the thresholds were trained heuristically to optimize sensitivity/specificity only in a single patient. In the current work, visual inspection of iEEG revealed that some of the events that fulfilled the Stockwell conditions were discarded in the first stage because they did not reach the threshold of mean  +  3  ×  SD. In addition, the same authors concluded that these three times of SD could be too restrictive and high. Table 5 shows the changes in precision, recall and F1 score values produced when the iEEG threshold was lowered to 2.5 and two times of SD. For this step, MEG-VS threshold values stayed invariant. An improvement in precision and F1 score values were produced when the iEEG threshold was lowered to 2.5  ×  SD, and only slightly lower recall values were observed. When the threshold was further decreased to 2  ×  SD, precision increased slightly but recall decreased more, producing a lower F1 score. These results suggested that iEEG threshold was fixed to be too restrictive for the data presented in this study and a slightly lower value would have been more appropriate.

The average number of detected ripples inside the ROI was 31.8  ±  13.2 and this value decreased to 13.4  ±  5.0 when computed outside the ROI (mean  ±  SD). The number of VS_WS involved in each ripple event was also computed. An average number of 18.4  ±  4.0 VSs were involved in each ripple event inside the ROI, whereas this value was slightly lower outside the ROI which showed an average of 12.6  ±  3.8 VS. The differences (inside and outside the ROI) between the number of detected ripples and the average number of VSs were statically significant (paired t-tests, significance set to 0.05) with probability values of 0.002 and 0.02, respectively. Table 4 shows the number of detected ripples outside the ROI and the mean number of VSs involved in each ripple for each patient.

The localization maps inside and outside the ROI for patient 1 are shown in figure 7. These maps measure the number of times that a specific virtual sensor was associated to a ripple. Figure 8 shows the grid inside and outside in slices (of the z-plane) for the nine subjects. These figures suggest that specific VSs appeared clustered, close to each other and involved in a high number of ripples inside the ROI, showing a more focalized localization map than outside the ROI. In order to quantify this effect, the percentage of detected ripples for the most active region inside and outside the ROI is shown in table 4. To compute this score, the VS showing the highest number of ripples and its six closest neighbors were selected for each subject. The ratio of the average number of ripples detected in seven clustered VSs to the number of detected ripples was then computed. As expected, the percentage of ripples involved in the most active region inside the ROI was significantly higher (p-value  <  0.001) than outside the ROI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are only a few studies that identify interictal epileptic HFOs using MEG in the time domain. In [23], a semi-automatic approach was proposed to detect FO (40–160 Hz) using MEG, where the proposed algorithm was adjusted to have high sensitivity, discarding FPs by human intervention (two experts). Visual detection of HFOs in the temporal domain is a time-consuming and highly subjective procedure. The current work represents the first study proposing an automatic approach that takes into account the time-frequency characteristics of the high-frequency events to detect ripples, while discarding other high-frequency noise events. Moreover, von Ellenrieder et al performed the detection of ripples at the sensor level and consequently they obtained low detection rates of ripples due to the low SNR of the signals at the scalp level. Thanks to beamformer-based spatial filtering, a high amount of high-frequency noise reduction was achieved, allowing to discern ripples that could not be identified in scalp MEG (figure 6). The remaining noisy events were discarded using a two-step algorithm based on thresholds whose values were selected to maximize the correspondence with ripple events in iEEG, the current gold standard. However, beamforming algorithms can be computationally consuming and a hybrid approach combining algorithms that detect some events at the sensor level (such as the one proposed by von Ellenrieder et al [23]) and algorithms that are focused in the source level could provide a good compromise between computation time and detection accuracy in MEG signals.

Furthermore, the whole-head beamformer step allowed choosing the ROI without relying on epileptic spikes. The exact relationship between IEDs and HFOs remains unclear, and it is unknown whether HFOs are produced by the same generators as IEDs or not [43]. Interictal HFOs have proven to be reliable markers of the SOZ, and better than epileptic spikes [43, 44]. In [24], the authors defined an ROI based on the location of the spikes in the physical MEG sensors, but this ROI was inaccurate in two out of 12 patients. In the current study, a finer VS grid was used and fixed to the area where most high-frequency events took place, minimizing the possibility of missing ripples occurring inside the volume of the grid.

To the best of the authors' knowledge, this is the first study that compares the simultaneous detection in time of ripples in noninvasive signals (MEG) against iEEG, the current gold standard for the detection of HFOs. Zelmann et al [18] had previously analyzed the spatial distribution between HFOs in iEEG and scalp EEG. In their study, the poor spatial correlation found between the two modalities was justified by the low number of channels (19-channel EEG system), the distortion produced by the electrical conductivity of the surrounding tissues and the low SNR of the signals. The use of systems with a higher number of channels and less influenced by the conductive properties of the skull, such as MEG systems, and methods for improving the SNR of the signals, such as beamforming, provides better results in terms of localization. The automatic analysis was performed in iEEG and the values for two thresholds were selected for MEG-VS to achieve a maximum F1 score; in other words, the weighted average of precision and recall measures. Unfortunately, the exact location of the iEEG in the head or the MRI of the subjects was unavailable for this study. Furthermore, the direct comparison of the recorded oscillations between these two signal acquisition methods is still uncertain. On the one hand, iEEG is an invasive technique that records electrical activity inside the head yielding a high SNR and consequently can detect events that could remain hidden in MEG signals, even after beamforming. Although there is no previous study evaluating HFOs in MEG and iEEG simultaneously in time, there are studies comparing simultaneous IED detection and their values of recall (or sensitivity) in MEG-VS detection with respect to iEEG, which were comparable to the results of the current study, around 60–80% [45–47]. While a more exhaustive comparison between the detection of IEDs and HFOs with both techniques under the same conditions should be carried out, the results of the detection should be comparable because IEDs and ripples have similar SNR values at their respective frequency bands [11], and because signals are equally attenuated at all frequencies of interest [18]. On the other hand, MEG can often provide additional and useful information mainly because of its whole-head coverage, while iEEG recordings usually only provide a locally limited neurophysiological picture [13]. The number of depth electrodes for this study covered a small region of the brain. For this reason, there is still some uncertainty regarding the FP ripples because they could be either noise being incorrectly detected as HFOs or actual ripples not recorded by iEEG.

The automatic detection of HFOs in invasive records is still an open field of study [9, 10]. For all these reasons, using the iEEG as the gold standard implies the acceptance of the detection provided by this technique. Table 5 shows that changing the SD threshold in the iEEG automatic detection produced variations in the number of co-detected events, suggesting that the threshold proposed in [10] might be fairly restrictive depending on the characteristics of the signals. The abovementioned selection of values for both thresholds considered the highest values of F1 scores, but they could be selected differently depending on the purpose of the detection: the values of factor_SD and $\text{thre}{{\text{s}}_{{{E}_{\text{r}}}}}$ could be decreased to maximize recall (higher than 80%) and to detect the highest possible number of ripples or, conversely, they could be increased to maximize precision (around 95%) and ensure almost totally that the detected events were equally visible in iEEG and MEG-VS. Nonetheless, more efforts should be directed to provide a reliable detection of ripples in MEG-VS by observing other characteristics of the actual signals (such as HFO onsets or connectivity between VS) or taking into account post-surgical outcomes (which were not available in this database).

An automatic selection of the area that showed the highest number of EOIs was performed to circumscribe the detection of ripples to a volume of approximately 9 cm³. This area agreed with the zone targeted by clinicians (tables 1 and 2). The HFO detection algorithm was applied to the same volume in the contralateral hemisphere with the purpose of comparing the number of ripples and the pattern drawn by VS_WS involved (figures 6 and 7). In general, a group of VSs that were detected simultaneously in different ripples produced a focalized map inside the ROI. Although this information was still insufficient to provide a reliable localization of the area generating pathological HFOs because of the lack of precise clinical information, it could be inferred that the oscillations were mostly generated by a common focus and this provides a general idea of the area where the ripples appeared. It is important to remark that although the automatic detection algorithm uses an ROI to detect ripples, these can be mapped a posteriori throughout the head by using temporal information about the detected ripple segments.

The proposed algorithm also detected some ripples outside the ROI, but the number of events was significantly lower and involved fewer VS. There is still uncertainty between the differential characteristics and mechanisms of pathological and physiological ripples [48]. Whether all the detected events were actually physiological, pathological or noise cannot be demonstrated definitely, but the detected ripples outside the ROI presented lower values of energy, which is a characteristic of physiological HFOs [49]. Moreover, the pattern of these oscillations appeared to be more scattered throughout the volume, indicating that these oscillations were not generated by a common area or focus. This finding is in correspondence to observations with animal studies indicating that the areas generating physiological HFOs appeared more extended [50] than the areas generating pathological HFOs, which are smaller brain regions on the order of cubic millimeters [51]. The results from Zelmann et al [18] suggested that it is possible to measure the activity from small cortical regions at scalp level. This could be explained by the solid angled concept introduced by Gloor [52]. These measurements should be comparable with those obtained with subdural grid contacts, always taking into account that the spatial resolution of the technique should be high.

Beamformer VSs allowed the detection of high-frequency activity otherwise not observable by physical sensors. The main purpose of this study was to automatically detect ripple events in MEG and discern them from noise by analyzing time-frequency characteristics using iEEG recordings to adjust and validate the algorithm. However, due to the lack of spatial location information for the intracranial electrodes, only a full-scale spatial validation of the detected ripples could be performed.

The spatial resolution of the reconstructed VSs is limited and depends on the depth, position and orientation of the sources, the SNR of the raw data, and the geometry of the head [53]. For all of these reasons, beamformer-based methods are not suitable for clinical application yet for the localization of the epileptic focus through ripples and the replacement of invasive techniques. However, automatic detection using non-invasive techniques could be useful to guide the implantation of the intracranial electrodes of subdural grids as it has been demonstrated that HFOs delimitate the affected area better than IEDs [43].

More effort should be made to increase the SNR of the signals, and in this sense independent component analysis has proven to be effective for improving the spatial localization of beamformer spatial filters [33]. Ultra-low-noise EEG/MEG techniques [54] are promising in the noninvasive detection of HFOs. Moreover, using the actual MRI of each subject to compute the head model would produce better results in the localization of the epileptic sources. A complete analysis comparing MEG with intracranial recordings should also include the position of the subdural electrodes.

Furthermore, using MEG systems with more channels and higher frequency rates would provide a better spatial and temporal resolution. Several studies have found that fast ripples have shown to be more correlated with the SOZ and presumably more highly epileptogenic [43, 55–57] than ripples. Due to the sampling rate limitation, this study only analyzed ripples up to 120 Hz. However, the automatic detection algorithm does not depend on the sampling rate and can be used to detect HFOs of higher frequency content acquired at higher rates. This would allow evaluating the differences in localization of ripples and fast ripples and their relationship with the SOZ. Furthermore, fast ripples would likely be downsampled to appear at lower harmonics, being able to detect them, but not to classify them between ripples or fast-ripples. Moreover, the study of ripples by non-invasive methods can also help to lateralize the epileptic focus in secondary bilateral synchrony [58], delimitate the SOZ in patients where fast ripples are not visible [55] and, in general, provide a better delimitation of the affected area than IEDs [43]. Moreover, detection of this area using non-invasive techniques could provide relevant information to guide the implantation of intracranial electrodes and for surgical treatment [23]. It is important to highlight that computation time increases with frequency and spatial resolution. These algorithms should be implemented efficiently, thereby minimizing the computation time to a level suitable for use in clinical analysis.

Finally, improving the spatial and temporal resolution would also allow detecting HFOs and determining the SOZ in the whole head using automated methods and clustering analysis, leading to computing different features such as density, connectivity, peak frequency, power, amplitude [59] or spectral entropy [9]. The automatic detection algorithm presented in this study would also be useful to extend these studies, which are currently limited for iEEG, and would exploit the advantage of the whole-head coverage of MEG, allowing assessment of the neural mechanisms of epilepsy as well as reducing invasiveness in clinical patients. A comparison between the areas of the brain generating HFOs and IEDs in epilepsy and their relationship would also be an interesting field of study. Moreover, a careful analysis of the whole brain volume could allow us to understand the areas involved in HFO generation, the propagation of this activity and how it is related with epileptogenic networks.