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

Objective. In epilepsy, high-frequency oscillations (HFOs) are expressively linked to the seizure onset zone (SOZ). The detection of HFOs in the noninvasive signals from scalp electroencephalography (EEG) and magnetoencephalography (MEG) is still a challenging task. The aim of this study was to automate the detection of ripples in MEG signals by reducing the high-frequency noise using beamformer-based virtual sensors (VSs) and applying an automatic procedure for exploring the time-frequency content of the detected events. Approach. Two-hundred seconds of MEG signal and simultaneous iEEG were selected from nine patients with refractory epilepsy. A two-stage algorithm was implemented. Firstly, beamforming was applied to the whole head to delimitate the region of interest (ROI) within a coarse grid of MEG-VS. Secondly, a beamformer using a finer grid in the ROI was computed. The automatic detection of ripples was performed using the time-frequency response provided by the Stockwell transform. Performance was evaluated through comparisons with simultaneous iEEG signals. Main results. ROIs were located within the seizure-generating lobes in the nine subjects. Precision and sensitivity values were 79.18% and 68.88%, respectively, by considering iEEG-detected events as benchmarks. A higher number of ripples were detected inside the ROI compared to the same region in the contralateral lobe. Significance. The evaluation of interictal ripples using non-invasive techniques can help in the delimitation of the epileptogenic zone and guide placement of intracranial electrodes. This is the first study that automatically detects ripples in MEG in the time domain located within the clinically expected epileptic area taking into account the time-frequency characteristics of the events through the whole signal spectrum. The algorithm was tested against intracranial recordings, the current gold standard. Further studies should explore this approach to enable the localization of noninvasively recorded HFOs to help during pre-surgical planning and to reduce the need for invasive diagnostics.

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Epilepsy is a neurological disorder that affects about 1% of the world population [1]. The most common therapeutic treatment relies on antiepileptic drugs (AEDs), which aim to reduce or eliminate seizures. Approximately 20% to 40% of patients diagnosed with epilepsy become refractory to treatment [2]. For these patients, the most frequent therapeutic alternative is neurosurgical resection of the epileptogenic zone (EZ), defined as the cortex area that is necessary and sufficient for initiating epileptic discharges and whose removal is necessary for complete abolition of seizures [3].
In current practice, the EZ is approximated by the seizure onset zone (SOZ), defined as the area where clinical seizures originate. Delimitation of this area is often done through presurgical analysis with intracranial electroencephalography (iEEG) [4] during the pre-ictal and ictal periods. In several recent studies, high-frequency oscillations (HFOs) have been evaluated as a specific biomarker for epileptogenicity in iEEG [5][6][7]. These fast oscillations (FOs) appear in the frequency range of 80-500 Hz, classified into ripples (80-200 Hz) and fast ripples (200-500 Hz), and are defined as spontaneous patterns above the baseline, clearly distinguished from noise and with at least four oscillations [8]. They can be recorded during the interictal period, without having to record ictal activity and, thereby, reducing discomfort and risk for the patient as well as recording time. HFOs may appear in concordance with interictal epileptiform discharges (IEDs) but also in isolation. The relationship between them is still an open discussion, although HFOs seem to be more specific to delimitate the EZ than IEDs [5].
Visual marking of HFOs is highly time-consuming. For iEEG, several automatic algorithms have been proposed in the time domain and in the time-frequency domain [9][10][11][12]. Time-frequency domain detectors can considerably reduce false positives (FPs), but they are also computationally more demanding. The detection of HFOs in non-invasive signals such as scalp EEG and magnetoencephalography (MEG) is still a challenging task, mainly because of the low signal-to-noise ratio (SNR) in high-frequency bands [13]. MEG and EEG signals must be adequately preprocessed to reduce biological and environmental interference coming from various sources [14]. Recent studies have shown that FO related to epilepsy can be identified in scalp EEG [15][16][17][18]; however, depending on the number of channels and due to the dist ortion of the signals caused by the skull, this technique may attenuate high frequency activity and present poor spatial resolution [18].
As it is much less influenced by the electrical conductivity of the surrounding tissues on the cerebral cortex [19], MEG can provide a more accurate non-invasive biomarker for epileptogenicity [13]. Various MEG studies have investigated high frequency bands [20][21][22]. There are few recent studies that have successfully detected high frequency activity in the time domain. In a study developed by von Ellenrieder et al [23], automatic detection of FO  previously developed for scalp EEG [11] was followed by a visual discarding of FPs performed in scalp MEG and then compared with the epileptogenic region defined by two specialists. However, presumably due to the low SNR of the FOs, the detection produced unsuccessful source localizations in some patients. van Klink et al introduced an interesting approach to detecting HFOs by improving the SNR using beamformer localization [24], one of the most used approaches for detecting MEG sources [25]. The beamformer reconstructs neuronal activity in a particular position as a weighted contribution from different sensors acting as a spatial filter capable of reducing and attenuating noise from distant sources [13,26,27]. Time-domain signals can be reconstructed at specific locations, usually known as virtual sensors (VSs), that can reveal information that is not discernable in physical sensors [24]. The detection of HFOs was done visually and the determination of the region of interest (ROI) was based on the localization of IEDs. This region is commonly known as the irritative zone [4]; it is usually more extensive than the theoretical EZ and it may even not overlap it.
In this study, an automatic method to detect ripples (80-120 Hz) in MEG is proposed. To remove high frequency noise, a beamformer analysis similar to [24] was performed in a two-step approach. VSs were used for the identification of epileptic HFOs that could not be assessed at the MEG sensor level. A simpler version of this approach was developed for a preliminary study [28] where HFOs were visually observed in the ROI and compared with simultaneous iEEG recordings. In this study, a comparison between the automated detection of ripples in MEG and iEEG was performed. To the authors' knowledge, this is the first automatic algorithm that can detect HFOs in MEG time-course signals.

Patients, signal acquisition settings and previous preprocessing
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][33][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).

Volume conduction model
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].

MEG automatic detection of ripples
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. 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): 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.

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.

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.

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.

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.

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:

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: 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.
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: 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 thres Er was defined to discard isolated events or those with low energy and, thereby, obtain better algorithm performance. Several values for thres Er were considered from 0% to 100% in steps of 1%.

Threshold selection: precision/recall curves
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 thres Er ) 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: 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 thres Er were selected using the median values for each of the leave-one-out functions.

Detection and localization of ripples inside and outside the ROI
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.

Detection of the ROI
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).

Threshold selection inside the ROI
To select the values of factor SD and thres Er 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 thres Er 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 thres Er : higher precision but lower recall with higher values.

Automatic detection
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-thanchance 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.  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.

Detection and localization of ripples inside and outside the ROI
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.

Main contribution of the study
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 timeconsuming 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 highfrequency 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 Table 4. Number of ripples, mean number of MEG-VS in each ripple for all the subjects and percentage of detected ripples inside the most active region. The evaluation outside the ROI was performed placing a grid of the same size and number of channels in the contralateral hemisphere where the ROI was placed.  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.

Threshold selection and performance analysis.
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][46][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 thres Er 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).

Evaluation of ripples inside and outside the ROI
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 3 . 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 differ ential 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.

Limitations of the study and further work
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 noninvasive 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][56][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.

Conclusion
HFOs are promising biomarkers for epilepsy. Several studies detect them automatically in intracranial recordings, but there are few studies evaluating these oscillations in noninvasive signals. To the authors' knowledge, this is the first study dealing with full automatic detection of ripples in MEG in the time domain taking into account the time-frequency characteristics of the events through the whole signal spectrum. Moreover, these ripples corresponded to the clinically expected epileptogenic area. The detection was compared with intracranial EEG recordings, demonstrating that the identification of ripples in MEG is feasible. The noninvasive study of HFOs is an interesting field of study that would allow evaluating HFO mechanisms in larger groups of epileptic patients and also in healthy subjects to discern the differences between pathologic and physiologic FO.