Mapping Interictal activity in epilepsy using a hidden Markov model: A magnetoencephalography study

Abstract Epilepsy is a highly heterogeneous neurological disorder with variable etiology, manifestation, and response to treatment. It is imperative that new models of epileptiform brain activity account for this variability, to identify individual needs and allow clinicians to curate personalized care. Here, we use a hidden Markov model (HMM) to create a unique statistical model of interictal brain activity for 10 pediatric patients. We use magnetoencephalography (MEG) data acquired as part of standard clinical care for patients at the Children's Hospital of Philadelphia. These data are routinely analyzed using excess kurtosis mapping (EKM); however, as cases become more complex (extreme multifocal and/or polymorphic activity), they become harder to interpret with EKM. We assessed the performance of the HMM against EKM for three patient groups, with increasingly complicated presentation. The difference in localization of epileptogenic foci for the two methods was 7 ± 2 mm (mean ± SD over all 10 patients); and 94% ± 13% of EKM temporal markers were matched by an HMM state visit. The HMM localizes epileptogenic areas (in agreement with EKM) and provides additional information about the relationship between those areas. A key advantage over current methods is that the HMM is a data‐driven model, so the output is tuned to each individual. Finally, the model output is intuitive, allowing a user (clinician) to review the result and manually select the HMM epileptiform state, offering multiple advantages over previous methods and allowing for broader implementation of MEG epileptiform analysis in surgical decision‐making for patients with intractable epilepsy.

offering multiple advantages over previous methods and allowing for broader implementation of MEG epileptiform analysis in surgical decision-making for patients with intractable epilepsy.

K E Y W O R D S
epilepsy, hidden Markov model, interictal activity, magnetoencephalography

| INTRODUCTION
Epilepsy is a neurological disorder affecting $50 million people worldwide (Beghi et al., 2019). It is characterized clinically by the occurrence of seizures which are generated by abnormal electrical cellular discharges in the brain. Epilepsy falls broadly into two categories: focal or generalized: in patients with focal epilepsy, there can be a single focus or multiple discreet foci that generate seizures, whereas for generalized epilepsy, seizures originate from diffuse areas throughout the cortex and/or deep brain structures such as the thalamus or hypothalamus. For many patients, antiepileptic medications can control both the severity and frequency of seizures. However, $30% of patients do not respond completely to medications (Mohan et al., 2018) and for these patients, surgical resection of affected brain tissue could be a viable treatment.
However, this necessitates extensive presurgical planning to accurately locate the affected brain area(s) prior to resection.
Current evaluation of pharmaco-resistant epilepsy is accomplished using electroencephalography (EEG; which measures electrical activity in the brain via assessment of electrical potentials at the scalp) alongside clinical factors and structural assessment using magnetic resonance imaging (MRI). This may be augmented by the nuclear medicine techniques of positron emission tomography and especially ictal (during seizure) single photon emission computed tomography. If candidate epileptogenic locations are identified via imaging, a more invasive intracranial EEG (iEEG) may be initiated where electrodes are located on the brain surface or within the grey matter. This allows electrophysiological assessment with optimal sensitivity and spatial accuracy, prior to resective surgery. Despite this extensive surgical planning, fewer than 50% of patients are seizure free 5 years postsurgery, with this number dropping to 38% at 10 years postsurgery (Mohan et al., 2018). It is therefore clear that a greater understanding of this heterogeneous disease, as well as improvements in clinical evaluation are required to improve patient outcome.
Magnetoencephalography (MEG) measures the magnetic fields induced by neuronal current flow. Unlike the electric potentials measured by EEG, magnetic fields are relatively unaffected by the high resistivity of the skull, resulting in less spatial distortion of the MEG compared with the EEG signal, and thus improved resolution. MEG is used in a growing number of clinical settings, particularly in epilepsy.
Not only does MEG provide additional information about the location of the epileptogenic zone (Agirre-Arrizubieta et al., 2014;Gofshteyn et al., 2019;Murakami et al., 2016;Nissen et al., 2016;Stefan et al., 2011), it can also be used to distinguish epileptogenic regions from eloquent cortex (Kim et al., 2013), and could be useful in mapping nonlesional or MRI negative focal epilepsies where there is no clear structural abnormality. Most importantly, a recent study by Rampp et al., (2019) showed, in 1000 patients, that presurgical MEG increases the chances of a patient achieving seizure freedom postsurgery when MEG localizations are resected.
MEG therefore has significant promise for assessment of patients with epilepsy, and the recent introduction of new technologies to capture the neuromagnetic field offer even higher spatial resolution, better sensitivity, and improved practicality, at lower cost (Boto et al., 2018(Boto et al., , 2021Brookes et al., 2022;Hill et al., 2020). This means that MEG could become even more established as the technique of choice for epilepsy evaluation. However, the detection of epileptogenic activity in MEG data remains a significant challenge. Recording data during a seizure is difficult due to uncontrolled patient movement, and for this reason, most MEG recordings are limited to interictal (between seizure) assessment (although see e.g., Tang et al., 2003).
Interictal events-sharply contoured atypical signals (known as epileptiform activity, i.e., spikes, sharps, etc.) are observable in resting MEG (and EEG) data and are generally assumed to originate from seizure onset zones, meaning that spatially mapping their origin offers useful information on the location of epileptogenic cortex. However, detecting interictal discharges can be challenging for two reasons. First, they are sporadic and unpredictable and, in some patients, rare. Thus, capturing interictal epileptiform activity can sometimes be a challenge without lengthy recording sessions (which should ideally include natural sleep because interictal discharges are often enhanced or only present in sleep). Even when they do occur, it can take significant time for a neurophysiologist to identify, mark, and categorize them (manually) in a MEG recording due to the high channel density relative to clinical EEG (i.e., >250 channels compared with <30).
Second, the temporal morphology of epileptogenic activity can vary markedly between patients; some produce spikes with/or without slow wave activity, which varies dramatically in amplitude. Other patients generate polymorphic bursts of "sharp wave" activity characterized by high temporal frequency signatures. In some individuals, the pattern of epileptiform activity repeats (like a template), in others it differs on each occurrence. This makes automatic detection algorithms challenging to design.
There are two commonly used analysis methods for localization of interictal activity, equivalent current dipole (ECD) fitting and excess kurtosis mapping (EKM). In ECD, interictal spikes are inspected visually (Bagi c et al., 2011). Once identified at the sensor level, a current dipole model is used to approximate the measured magnetic field just prior to, or at the peak of, the spike; by letting the origin of the modeled dipole vary spatially, and then selecting the point at which the model best fits the measured field, it becomes possible to localize the brain region generating the spike (Ebersole, 1997;Wheless et al., 1999). This technique works reasonably well in cases where high amplitude spikes are observed in isolation but is often not useful in cases where interictal activity includes polymorphic bursts. In addition, multifocal epilepsies are a challenge since ECD requires a priori estimation of the number of active regions.
In contrast, EKM (Gaetz et al., 2015;Robinson et al., 2004;Schwartz et al., 2010) is an automatic method which assumes only that the epileptiform signals of interest are sharply contoured (relative to the typically rhythmic background MEG signals). Kurtosis is a measure of the shape of a statistical distribution; in cases of abnormal activity (e.g., if a dataset has large spikes), it's statistical distribution includes a large "tail" and thus begins to look non-Gaussian; this is quantified by increased kurtosis (also known as the fourth moment of the distribution). By application of a kurtosis algorithm to MEG signals extracted from multiple brain locations, it becomes possible to localize areas generating abnormal activity. EKM does not require a priori estimation of the number of epileptogenic regions; further, it is not limited to spikes, but can be used to assess any atypical activity, provided it has high kurtosis. However, EKM also has limitations; it has low sensitivity to low-amplitude polymorphic activity. Also, counterintuitively, in cases with rapidly occurring high-amplitude spikes, excess kurtosis has diminished sensitivity, because the kurtotic signals are so common they begin to represent the mean. For these reasons, neither ECD nor EKM is a perfect solution to analysis of MEG data in epilepsy, and other methods, which can accurately and automatically identify epileptiform activity, map its spatial origins, and (in multifocal epilepsy) characterize relationships between regions, would be useful.
Hidden Markov modeling (HMM) has gained traction in recent years as a method to elucidate complex neural dynamics in MEG data (Higgins et al., 2021;Quinn et al., 2018;Vidaurre et al., 2018). The method works by detection of repeated patterns of activity (known as states) in MEG data; patterns can be characterized based on a number of features, including amplitude, channel covariance, and spectral properties. In this article, we aimed to test the hypothesis that "an HMM could be used to identify, in time and space, epileptiform activity in agreement with current state of the art methods". We further aimed to show that, "in at least some patients, our method could offer more information than the established EKM technique". In what follows, we describe our data processing pipeline, and its application to MEG data acquired in multiple different patient groups, ranging from "simple" cases with focal spike-and-wave epileptiform activity for a single locus, to more complex cases with multifocal epileptiform activity exhibiting polymorphic bursts from multiple loci.

| Patient identification and data collection
This study was determined by the IRB to have exempt ethics status as it constituted secondary analysis of data captured for clinical purposes under the NIH common rule (January 2019). All research subjects were scanned as part of clinical care at the Children's Hospital of Philadelphia. Data from 10 pediatric patients undergoing presurgical epilepsy evaluation between the ages of 11 months and 17 years (median 8.18 years, see Table 1), were utilized in the study. There were two female and eight male patients. In total, 9 of the 10 had focal epilepsy with or without impaired awareness, and one (Patient 6) had epilepsy with combined focal and generalized features. Six patients had right focal epileptiform discharges, whereas 4 had left focal epileptiform discharges, as identified by clinical features and multimodal imaging methods including MEG analyzed with EKM.
Most focal discharges were in the frontal or temporal regions with one patient (Patient 8) with central localization and one patient with both frontal and posterior temporal discharges (Patient 6). When known, etiology was primarily confirmed or suspected structural abnormality (such as focal cortical dysplasia), with one patient (Patient 1) with confirmed genetic etiology and two patients with, as yet unknown etiology.
Multiple 2-min recordings were acquired using a CTF 275-channel MEG system operating in third-order synthetic gradiometer configuration. By collecting data in relatively short 2-min runs, the maximum head motion within each dataset was minimized. The number of acquired datasets per patient ranged from 15 to 29; mean = 17.8. Data were acquired at a sample rate of 1200 Hz. In all cases, patients were scanned supine. Of the 10 patients, 7 were scanned whilst sedated using general anesthesia.
Prior to data acquisition, three head position indicator coils were attached to fiducial points on the head. During recording, these coils were energized (at nonphysiological AC frequencies) to allow continuous localization of their position relative to the MEG sensor array. All MEG scans were followed by an anatomical MRI during which MRI contrast markers were placed at the same fiducial points on the head.
Coregistration of the MEG and MRI fiducial locations thus enabled complete spatial mapping of the MEG array relative to individual brain anatomy. This coregistration, in turn, allowed generation of functional images showing the cortical origins of epileptiform activity.
Given that the purpose of this work is to act as a feasibility study to test whether the HMM could be a useful tool in the identification and localization of epileptiform activity, patients were selected retrospectively to fall into three groups of increasing localization difficulty (see Table 1); four had focal epilepsy generating interictal spike and wave activity; two had focal epilepsy with only polymorphic bursts, and the final four had multifocal epilepsy with polymorphic bursts. To assess the efficacy of our HMM MEG method, results from iEEG and/or surgical outcomes are provided for patients where these data exist (Patients 2, 5, and 8). A comparison with scalp EEG is also given in the single most likely electrode), with the electrode positions of observed EEG discharges. Where our HMM MEG method was subsequently found to be concordant with EKM, it would therefore also be concordant with scalp EEG.

| Data processing
Following MEG collection, data were visually inspected and any 2-min runs containing obvious interference or segments where patient motion exceeded 1 cm were removed from further analysis. Since some patients have infrequent interictal activity, all data were checked by a MEG expert (William Gaetz) to ensure they contained epileptiform activity and runs which did not were subsequently discarded. This left an average of 11 runs per subject. Our HMM-based method of mapping epileptiform activity comprised a two-step process (see Figure 1): In step one, the HMM was applied to channellevel MEG data to identify time periods in which interictal epileptiform activity occurred. In step two, a beamformer was used to localize the brain regions generating the activity identified by the HMM.
These two steps are described in detail below.

| Hidden Markov modeling
To find spatiotemporal patterns corresponding to epileptogenic activity, we applied a multivariate, five-state, time delay embedded HMM, in channel space. An HMM assumes that a series of recurring mutually exclusive "hidden" states govern the MEG data, such that every point in time is associated with one of the states. The sequence is assumed to be Markovian (i.e., the state active at a time point, t, only depends on that active at time point t À 1). An observation model links the HMM state to the observed values in the MEG data.
The HMM has been described extensively in previous papers (Baker et al., 2014;Vidaurre et al., 2018) and a complete mathematical description will not be repeated here. Briefly, in its simplest form, an HMM would describe each state using a multivariate distribution; that is, a mean (for all channels) and covariance (across channels). The five distributions that best described the data would be derived, and the probability of each data-point belonging to a specific state would be calculated. The number of states is defined a priori and model inference would learn the sequence of states, from the observed data.
Here, we employed a more complex model which also allowed time-delay embedding , adding information on autocovariance (defined over a specified time window [duration 73 ms]). These state autocovariance patterns contain the spectral content of the signal, consequentially using our HMM, a single state is defined based upon signal variance, covariance across channels, and spectral content. This model had the potential to characterize both spike and wave activity and polymorphic bursts; in the former case the (typically high) amplitude of a spike, with a full width at half maximum of $70 ms, coupled with its distinct spectral content would characterize the state and differentiate it from ongoing "normal" activity. In the latter case, since polymorphic bursts are associated with "sharps" (high-frequency activity) we again reasoned that distinct spectral content would define the state.
Prior to application of the HMM, the data used for the model inference (comprising 266-channel sensor space MEG data) were F I G U R E 1 Schematic representation of the hidden Markov model (HMM)-based process to identify epileptogenic activity. A multivariate time delay embedded HMM was used to identify five states, each state characterized by its mean, covariance (across channels), and spectral content. A beamformer was used, along with temporal state allocations, to generate images of state activity across the cortex and a time course of activity from the peak voxel. This allowed us to identify an epileptiform state and a map of epileptogenic cortex. LCMV, linearly constrained minimum variance; MEG, magnetoencephalography bandpass filtered between 20 and 70 Hz (to match the standard EKM pipeline used by the Children's Hospital of Philadelphia-see also below), notch filtered at 60 Hz (to remove mains frequency artifacts), temporally down sampled to 150 Hz, time-embedded using 73 ms lags and a principal component analysis was used to reduce the data to 50 components (this allowed for faster model inference and helped to avoid overfitting). A total of fifty principal components were enough to capture 87% ± 3% of the data variance (mean and SD over subjects and runs). Note that each 2-min clinical run was considered separately, just as the EKM data were. The model inference itself was undertaken using a variational Bayesian method which seeks to minimize the free energy of the system. We computed five-states; for each state, in addition to an observation model, we obtained a time course of the probability that the state is active. These time courses were thresholded at two thirds, so the state was defined as "active" when the probability exceeded two thirds (Seedat et al., 2020). From the time courses, we also obtained a state-transition-matrix-a 5 Â 5 matrix of probabilities defining the temporal relationship between states (i.e., element 2,1 in the matrix would represent the probability that State 1 followed State 2). This approach might offer useful information in cases with multifocal epileptiform activity where one source consistently precedes the other.
To assess the impact of the choice of model parameters used for the HMM inference, the model inference was computed multiple times for various lag durations (ranging from 33 to 207 ms) and for various numbers of states (ranging from 3 to 9). Please see the Supplementary Information S1 for further details.

| Beamforming
Following application of the HMM (in sensor space), a linearly constrained minimum variance beamformer (Robinson and Vrba, 1999) was used to localize the spatial signature of each state in the brain.
The brain was divided into a regular 4 mm grid of voxels, and each voxel time course was defined as a weighted sum of sensor measures.
The beamformer weights were defined using a data covariance matrix calculated in the 20-70 Hz frequency band, and a time window spanning the entire recording. To maximize spatial resolution, no regularization was applied (Brookes et al 2008 every voxel time course in the brain, to highlight the brain regions which elicit changes in variance when the state switches on. This produced a spatial map of state activity. In addition to the spatial maps, we also used the beamformer to derive a time course of electrophysiological activity at the peak location of the spatial map. Beamformer weights were defined as above but using covariance calculated in the 1-150 Hz band, and a single time course showing 1-150 Hz activity was extracted. Having derived a spatial map and time course of activity for each of the five states, these were visually inspected by a single MEGepilepsy expert (William Gaetz). Those states whose time courses showed epileptiform activity when the state was active were identified, and the spatial localization was noted for each run. These were termed the "epileptiform state(s)."

| Comparison to existing methods
We compared the results of our HMM, to the more established EKM technique. We selected EKM for this comparison because of its advantages over ECD (Hall et al., 2018) and its use in large pediatric cohorts (Gofshteyn et al., 2019). To ensure that a standard EKM pipeline was followed, we used commercial software developed by CTF (Coquitlam, BC, Canada) known as SAM (g2), and the established pipeline used clinically by the epilepsy team at the Children's Hospital of Philadelphia (Schwartz et al., 2010).
Prior to the application of EKM, all data were filtered 20-70 Hz.
In the SAM(g2) implementation, the brain was divided into a regular grid of 5 mm voxels and a scalar beamformer (equivalent to that    The patient exhibited abundant, large amplitude spikes, with some slow wave activity. Spatial correspondence between EKM and HMM was 6 ± 2 mm (mean and standard deviation over eight 2-min runs) and 102 ± 8% of EKM-identified epileptic events were matched in time by the occurrence of the epileptiform state. Note that more than one state visit to a single kurtosis marker will occasionally yield values >100%. Across all datasets, the HMM epileptiform state was active for 4% ± 2% of the time.
iEEG and postsurgical outcomes were available for this patient, with grids over the right lateral frontal, right orbitofrontal, and right interhemispheric fissure, as well as with a depth electrode in the right lateral frontal region. Interictal activity was seen in the right lateral frontal grid and right frontal depth electrode but was also seen broadly in the other two grids. Seizures arose mainly from the right lateral frontal grid and depth electrode but also in the other regions almost simultaneously. The iEEG and MEG are therefore concordant, especially when considering the broad onset in the intracranial study and the MEG report. Although motor strip was involved, the patient had a right frontal lobectomy sparing motor strip which resulted in diagnosis of focal cortical dysplasia (FCD) type 2B. The patient was then seizure free for 2 years with recurrence in the motor region but has been seizure free on adjusted medication for 3 years since then. between the HMM and EKM. In these cases, spatial correspondence was 4 ± 1, 9 ± 5, and 3.9 ± 0.3 mm, and the temporal correspondence was 97% ± 4%, 100% ± 1%, and 96% ± 5%. In general, these results support our hypothesis that the HMM performs similarly to EKM in enabling the identification of epileptiform activity in time and space, at least in patients with focal spike and wave epileptiform activity.
Spatiotemporal correspondence was high in all cases and data showed good correspondence across many runs, for each subject.

| Case 2: Focal epilepsy with polymorphic activity
Case 2 ( Figure 4) shows MEG data acquired in a patient with focal epilepsy, but without typical spikes in the MEG trace (Figure 4b). The patient's resting MEG data exhibited abundant polymorphic bursts of sharply contoured epileptiform activity, which (unlike spikes) change their temporal morphology on each occurrence. Such data are not amenable to conventional ECD source analysis; however, here we see that both the HMM and EKM generate a focal localization of the epileptogenic zone with excellent spatial agreement between the two methods. On average across fifteen 2-min recordings, the spatial discrepancy between the HMM and EKM peak location was 5 ± 2 mm.
In addition, the temporal coincidence of the epileptiform HMM state and the EKM-derived markers was 106% ± 7%. The percentage of time when the HMM epileptiform state was active was 6% ± 3% (average and SD over all runs).
iEEG for this patient was stereo-EEG. Interictal activity arose mainly from the mid superior temporal gyrus (STG) electrode, as well F I G U R E 3 Epilepsy Case 1magnetoencephalography (MEG) localizations compared with resected area for patient 2 (focal spike wave). Panel (a) The MEG localizations for the excess kurtosis mapping method (red markers) and the hidden Markov model method (yellow markers). Panel (b) shows the resected area in the postoperative magnetic resonance imaging (MRI). Note that the postoperative MRI has different resolutions in the sagittal, coronal, and axial planes as seizures, but there were rare interictal discharges in the angular gyrus as well. The patient had a focal resection of the posterior STG which resulted in pathology of FCD 2B. She has been seizure free since surgery (9 months) and has been able to reduce medication. The surgery did result in a conductive aphasia, but this has since improved.
MEG was concordant in the post-STG but also had some involvement A second case of focal epilepsy with polymorphic bursts is given in Figure S6; results again were similar with a spatial correspondence of 5 ± 3 mm and temporal correspondence of 113% ± 14%.

| Case 3: Multifocal epilepsy with polymorphic activity
The above results show relatively straightforward cases of focal epilepsy, where abnormal epileptiform activity arises from a single location in the brain. However, in more complex cases, abnormal activity can occur from more than one location (often simultaneously), and in such cases, the challenge becomes determining where these multiple regions are, how they are related, and if possible, which region serves as the driver of an epileptiform network causing other areas to exhibit epileptiform activity.
Case 3 is a patient with multifocal epilepsy; in total 15 datasets were acquired in this individual, results from a single representative run are shown in Figure 6. This patient exhibited EKM and HMM peaks in right periorolandic areas. For 11/15 runs, there was just one epileptiform state describing the activity from these regions, in the other four runs two epileptiform states were identified. Where more than one HMM state was identified, there was no obvious temporal relationship between them.
In this patient, iEEG was grids and strips over the right superior active was 3% ± 1% and 4% ± 1% for the temporal and frontal epileptiform states, respectively.
We also compared the results shown in Figure 8 with the output of our EKM algorithm. EKM generated a single map which also had peaks in frontal and temporal lobe. On average (across 15 runs) the mean spatial correspondence for the frontal lobe peak was 8 ± 3 mm and the equivalent distance for the temporal lobe peak was 7 ± 4 mm, once again implying spatial correspondence between the HMM and EKM. This is impressive given the challenges posed by such a complicated case to each of these methods. It is worth noting that although the EKM method places temporal markers in the data to help epileptologists assess whether a spike in one part of the brain precedes a spike in another area, the HMM provides additional information about the temporal relationship between the two brain locations using the state transition probabilities. Furthermore, the HMM uses all of the state data to estimate transitions, something which becomes particularly important in cases without clear spikes. No invasive assessment data or surgical outcomes are available for this patient.
Two further cases of multifocal epilepsy are presented in Figures S7 and S8. In both cases, results are similar to those shown in Figures 6 and 8.
F I G U R E 5 Epilepsy Case 2magnetoencephalography (MEG) localizations compared with resected area for Patient 5 (focal polymorphic). Panel (a) shows the MEG localizations (peak virtual electrode locations) for the excess kurtosis mapping method (red markers) and the hidden Markov model method (yellow markers). Note that markers may be overlaid so that red markers are hidden under yellow ones. Panel (b) shows the resected area in the postoperative magnetic resonance imaging

| Group results
For each of the three patient groups, the average spatial and temporal correspondence between the two methods was found over subjects and this is shown in Table 2. There was good agreement between the two methods for all three groups. This is especially encouraging given the complexity of the multifocal patients, with less than a centimeter discrepancy between the peak locations. There was more variation in the temporal coincidence metric which is likely to be because there were fewer EKM markers in the multifocal polymorphic data because the amplitude of the virtual electrode data rarely exceeded the threshold needed for marker placement. This meant that if the HMM missed a single marker, it resulted in a much-reduced percentage coincidence.

| DISCUSSION
Epilepsy is a debilitating disorder in which both symptoms and treatments differ markedly across patients. In some cases where pharmacological intervention fails to control seizures, patients become candidates for surgery in which affected brain regions are resected.
For many patients, such intervention offers seizure freedom (and thus a marked improvement in quality of life). However, success greatly depends on presurgical planning to accurately identify the epileptogenic region(s) and the current clinical pathway (involving multidisciplinary team input with interictal and ictal EEG ± iEEG, and MRI) is not always successful in identifying candidate regions for resection (or indeed implantation of the iEEG). Consequently, improvements to this pathway could enable more patients to become eligible for surgery and might offer improved outcomes for those who do have surgery.
MEG has a significant and established role in epilepsy, offering high-precision mapping of epileptogenic and eloquent cortex (Kim et al., 2013;Rampp et al., 2019;Schwartz et al., 2010). However, the current methods for analysis of MEG data are limited. ECD-still the most widely used technique-is unsuited to multifocal epilepsies, and cases where temporal morphology of interictal events fails to include isolated high amplitude spikes. Some of the limitations of ECD are lifted via the use of kurtosis-based techniques; however, these too offer limited sensitivity; they only respond to activity that is characterized by non-Gaussianity, and in multifocal cases, EKM provides limited information about the relationship between brain regions. This F I G U R E 6 Epilepsy Case 3-Patient 8; multifocal, polymorphic. There was good agreement between methods for the spatial localization (see the maps for hidden Markov model [HMM; a, left] and excess kurtosis mapping [EKM; a, right]) with an average distance between HMM and EKM peaks of 11 ± 2 mm. The virtual depth electrode time courses for each method are shown in (b) with HMM state visits highlighted in red and EKM markers as dashed lines. The temporal match was also good with 88% ± 15% of EKM markers matched by an HMM state visit. The HMM state was active for 7% ± 3% of the total time means that identified brain areas can be erroneously discarded due to a lack of evidence that they are involved in an epileptic network. Consequently, there is room for more adaptable processing pipelines to analyze MEG data in epilepsy patients.
In this work, we have demonstrated that an HMM is a promising technique to localize epileptiform activity in both time and space. The HMM has been used increasingly in the analysis of MEG data in recent years, to elucidate the complex spatiotemporal dynamics of brain networks (Baker et al., 2014;Seedat et al., 2020;Vidaurre et al., 2018;Woolrich et al., 2013). Here we show that, by seeking short segments of data with unique spatiospectral characteristics, the HMM allows automatic detection of epileptiform activity. Importantly, the HMM does not rely on a single data characteristic-for example, high amplitude spikes, or a non-Gaussian statistical distribution. In addition, activity does not have to repeat in time (like a template).
Rather, the HMM employs features of the data (e.g., sharp contouring) that can vary on every occurrence of an epileptiform state. Consequently, as we have demonstrated, it can adapt to the heterogeneity of the disorder; the same algorithm will find spike and waves, sharps, or polymorphic bursts with equal efficiency. The flexible, data-driven approach of the HMM lends itself naturally to personalized medicine, where the output of the model will be specific to each patient.
We compared the output of our HMM to the more established EKM technique (which is embedded in a number of clinical pathways for epilepsy across multiple centers). Spatial agreement was within 8 ± 3 mm (mean and SD across all subjects and runs in all three patient groups), showing that epileptic foci could be identified using the HMM with high spatial acuity. Temporal agreement was also good, with 94% ± 13% (mean over all subjects in all patient groups) of epileptiform events detected by the standard EKM pipeline matched by the occurrence of the epileptiform state, identified by the HMM.
Importantly, whilst the HMM looks for a specific spatiospectral pattern, the EKM pipeline is only able to find the timing of epileptic bursts via the use of a threshold: specifically, the EKM threshold is set by some arbitrary parameter (6) multiplied by the RMS amplitude.
Thus, epileptiform activity with lower amplitude is likely to be missed by this overly simplified technique. The HMM on the other hand relies on the occurrence of a specific multivariate spatiospectral pattern.
Consequently, it is less likely to miss epileptiform bursts.
The HMM has demonstrated advantages when looking at multifocal epilepsy cases. In such complex cases, ECD is limited since one must make an a priori decision on the number of active foci, which is always unknown. EKM is better, since it allows identification of multiple regions, and the manual comparison of kurtosis markers can identify whether one sharply contoured event from one location consistently precedes another. Similarly, the HMM can identify multiple epileptiform foci but provides more information regarding the temporal relationship between the identified regions. Specifically, where multiple foci appear within a single state, then those regions (by definition) are related temporally, that is, the regions are "active" together. This was seen in Patient 9 ( Figure S8) in which (in 6/13 cases) an epileptogenic focus close to an area of resected tissue in the frontal lobe appeared alongside a second peak in parietal cortex. The signals from the second peak were not epileptiform and thus had initially been dismissed; however, the fact that the activity generated increases in variance in synchrony with the peak close to the resec- The HMM method implemented in this study is designed to be a semiautomated spike-detection algorithm to aid clinicians in quickly and efficiently identifying epileptiform discharges in temporal data as well as localizing the regions of the brain in which they arise. We anticipate that in the clinic, there will be two manual processing steps required: the first is the visual inspection of data to remove those datasets with excessive head motion. This could be easily undertaken by somebody without training in reading clinical MEG data or identifying epileptiform discharges and could potentially be automated (at least in part). The second step will be to interpret the output of the HMM to select the epileptiform state(s)-this will likely require the assistance of a trained epileptologist. This method utilizes both the power of machine learning (with the HMM) and clinical expertise, but there is the potential to further automate the pipeline by employing additional machine learning techniques to learn features of epileptiform states for automatic state selection in the future.
There are a number of further extensions to the present methodology which should be considered. First, the current version of the model operates using sensor space data, and localization is achieved using a beamformer. It is theoretically possible to apply the beamformer first, and then the HMM, and this may offer advantages: For example, the beamformer is known to supress fields from non-brain sources, which would mean a source space HMM is less likely to be biased by interference. However, it is also significantly more computationally demanding since the HMM must be run on many thousands of (voxel-based) signals rather than a few hundred channels. This was found to be impractical with current computing power, and whilst brain parcellation (e.g., dividing the brain into a smaller number of anatomically [or functionally] meaningful regions) offers a means to reduce dimensionality, it also reduces spatial specificity which is of key importance in this application. Thus, we feel the method presented is the most practical currently, but the use of better computing may ultimately enable a source space HMM. In addition, although we have applied our method to MEG data, application to EEG data or concurrent EEG/MEG data would also be valuable. In the former case, the limitations surrounding the spatial resolution of EEG would mean low spatial specificity, but nevertheless, the HMM might offer a useful means to automatically identify time points which contain interictal events-this could enable a significant saving in time for epileptologists who review EEG data. For concurrent EEG/MEG, there is some evidence that the combination increases sensitivity (e.g., to radial sources) This may offer a significant advantage over MEG alone. It should also be noted that the data used in this study were from a single site and acquired using one type of MEG system (CTF)-future studies should test the reliability of the HMM optimized with these data at other sites.
The advent of wearable MEG systems, based on optically pumped magnetometer detection of the neuromagnetic field Hill et al., 2019) may also offer significant advantages.
OPMs measure magnetic fields much closer to the brain than SQUIDs (superconducting quantum interference devices), offering greater spatial precision and sensitivity. Likewise, bespoke helmets, which can be removed and replaced, with sensors going back to the same locations relative to brain anatomy each time, may offer the opportunity to record longer datasets. This would undoubtedly create a more accurate model which could help to overcome the variations in HMM output over multiple runs. Wearable MEG would also enable free movement whilst scanning; one of the key limitations of our study is that seven of the patients were under general anesthesia to stop them moving. Unfortunately, anesthesia (even lightly applied) is likely to reduce all neural activity including epileptiform activity, and so scan-

| CONCLUSION
The HMM offers an alternative to ECD and EKM for identifying and localizing epileptiform activity in the human brain. In 10 subjects, we found good spatial agreement between methods, with the HMM able to provide localization which matched, within a few mm, that from EKM. Moreover, we found that the HMM offers more information about epileptiform activity, particularly in multifocal cases. As the use

DATA AVAILABILITY STATEMENT
Data available on request due to privacy/ethical restrictions.