Systems and method for multimodal epileptogenic zone identification and cognitive function mapping

Abstract

A comprehensive system includes a multimodal integrated database storing population-level data associated with brain activity, and various models for generating artificial intelligence (AI) models for detecting biomarkers in the data, and a training the models for feature-based classification. The system also includes a module for receiving multimodal data associated with a patient's brain, as well as a module for aggregating the patient data and trained AI models to generate brain maps and treatment plans. The system may be iteratively refined with outcome data or input by a user via a user interface. The user interface may also communicate the treatment plan to a user, which may include resection of one or more target areas in the patient's brain or request for further data acquisition.

Description

SYSTEMS AND METHOD FOR MULTIMODAL EPILEPTOGENIC ZONE IDENTIFICATION AND COGNITIVE FUNCTION MAPPINGCross Reference to Related Applications

[0001] The present application is based on, claims priority to, and incorporates herein by reference in its entirety for all purses, US Provisional Application Serial No. 63 / 508,245 filed June 14, 2023.Statement of Government Support

[0002] This invention was made with government support under NS 128318 awarded by the National Institutes of Health. The government has certain rights in the invention.Background

[0003] An estimated 1.2 million Americans suffer from epilepsy, and globally the number of such patients is around 1% of the population. Approximately 30% of such patients - around 25 million people globally- suffer from epileptic seizures that do not respond to drugs. Treatment of such drug resistant patients remains a very challenging task. There is more hope for around 60% of such patients (or around 15M people globally) who suffer from focal seizures, also known as partial seizures, which originate in just one part of the brain, and then spread across to other regions causing observable seizures. Locating epileptogenic tissues (that cause seizures to nucleate) provides at least two standard treatment options: (i) electrodes are implanted in the brain to conduct stimulation (intervention) at or near the locations with epileptogenic tissues to control or stop seizures before they spread, or (ii) when the associated regions would not lead to severe cognitive impairment, those regions can be resected via surgery, potentially leading to post-surgery seizure freedom. Both procedures involve considerable risks. For example, electrodes might get infected or become isolated by localized inflammations, leading to a loss of intervention capability. Similarly, resection can unintentionally remove tissues that are cognitively important leading to severe disabilities.

[0004] Thus, there is a need for high precision localization of epileptogenic tissues in the brain and a determination of cognitive functionalities of those brain regions.Summary

[0005] The present disclosure provides systems and methods that overcome the aforementioned drawbacks by providing a system-wide reinforcement learning system to converge on a treatment protocol for different patient conditions. The multi-module system may include a (1) multimodal patient database, (2) an artificial intelligence (Al) modelgenerator for patient-specific brain region localization and an Al-based patient-specific brain region inference and treatment recommendation for an individual patient, and (3) a user interface for providing input data, configure analysis settings, and view model outputs.

[0006] In one aspect of the present disclosure, a system for identifying pathological features in medical data is presented. The system comprises a computing device comprising a memory storing instructions and multimodal datasets and a processor configured to access the memory to execute the instructions. The processor is caused to receive one or more datasets from one or more modalities associated with brain activity of a population, identify one or more biomarkers for each of the one or more datasets associated with seizure activity, extract one or more features from the one or more biomarkers, identify cognitive spatial biomarkers for each of the one or more datasets, and perform a deep learning (DL)-based classification of the one or more biomarkers based on the one or more features. The computing device further includes a user interface configured to output at least one of the one or more biomarkers or the cognitive spatial biomarkers.

[0007] In another aspect of the present disclosure, a system for identifying treatment target regions of a patient based on pathological features in acquired patient data is presented. The system comprises a computing device comprising a memory storing instructions and a processor configured to access the memory to execute the instructions. The processor is caused to receive one or more features extracted from one or more biomarkers associated with seizure activity identified from one or more datasets from one or more modalities associated with brain activity of a population, obtain input data of a patient’s brain comprising one or more modalities to identify the one or more features in the input data, generate a brain map of the patient based on the one or more features identified in the input data, and generate a treatment plan for one or more target regions in the patient’s brain based on the brain map. The computing device further includes a user interface configured to receive the input data and output at least one of the one or more biomarkers, the brain map, or the treatment plan.

[0008] In one aspect of the present disclosure, a method for identifying pathological biomarkers in medical data is presented. The method comprises receiving, by a system comprising a processor, one or more datasets from one or more modalities associated with brain activity of a population, identifying, by the processor, one or more biomarkers for each of the one or more datasets associated with seizure activity, extracting, by the processor, one or more features from the one or more biomarkers, identifying, by the processor, cognitive spatial biomarkers for each of the one or more datasets, and performing, by the processor, a deeplearning (DL)-based classification of the one or more biomarkers based on the one or more features. The method further comprises outputting, by a user interface, at least one of the one or more biomarkers or the cognitive spatial biomarkers.

[0009] In another aspect of the present disclosure, a method for identifying treatment target regions of a patient based on pathological features in acquired patient data is presented. The method comprises receiving, by a system comprising a processor, one or more features extracted from one or more biomarkers associated with seizure activity identified from one or more datasets from one or more modalities associated with brain activity of a population, obtaining, by the processor, input data of a patient’s brain comprising one or more modalities, identifying, by the processor, the one or more features in the input data, generating, by the processor, a brain map of the patient based on the one or more features identified in the input data, generating, by the processor, a treatment plan for one or more target regions in the patient’ s brain based on the brain map, and receiving input data, by a user interface. The method further comprises outputting, by the user interface, at least one of the one or more biomarkers, the brain map, or the treatment plan.

[0010] These aspects are nonlimiting. Other aspects and features of the systems and methods described herein will be provided below.Brief Description of the Drawings

[0011] The foregoing features of embodiments will be more readily understood by reference to the following detailed description, taken with reference to the accompanying drawings, in which:

[0012] FIG. 1 A is a flowchart of the framework, according to aspects of the present disclosure.

[0013] FIG. IB is a schematic of the model generation of the framework of FIG. 1 A, according to aspects of the present disclosure.

[0014] FIG. 1C is a schematic of the training plan generation of the framework of FIG. 1A, according to aspects of the present disclosure.

[0015] FIG. 2A is a flow chart of the method of generating and training Al models, according to aspects of the present disclosure.

[0016] FIG. 2B is a flow chart of the method of determining a treatment plan for a patient, according to aspects of the present disclosure.

[0017] FIG. 3A is a flow chart of an example implementation of an inference engine for processing biomarkers in EEG signals, according to aspects of the present disclosure.

[0018] FIG. 3B is a flow chart of an example implementation of an inference engine for processing biomarkers in EEG signals, according to aspects of the present disclosure.

[0019] FIG. 4A is an example implementation of a graphical user interface, according to aspects of the present disclosure.

[0020] FIG. 4B is a subset of options available in the graphical user interface of FIG. 4A, according to aspects of the present disclosure.

[0021] FIG. 5 is a comparison between Matlab (left) and Scipy (right) regarding filter construction. For the given parameters FP= 80 Hz, Fs= 500 Hz, TP= 0.5 dB, and Ts= 100 dB, which is specified in RIPPLELAB, Matlab did not precisely match the frequency response to the desired PassBand Ripple and StopBand Attenuation. In contrast, Scipy utilized in our implementation generates a more closely aligned frequency response. Owing to numerical overflow, we use a slightly adjusted StopBand Attenuation value of Ts = 93 dB in our PyHFO implementation (Fp:PassBand, Fs:StopBand, Tp:PassBand Ripple, Ts: StopBand Attenuation).

[0022] FIG. 6A is a plot of the time-domain augmentation consisting of two steps: 1) random time-domain shift, the window of selecting EEG tracing to generate feature (initially center at the middle ofthe HFO with ± 285 ms around it, purple region) could be randomly shifted ± 50 ms in the time domain (orange region). 2) the EEG signal was randomly flipped in the time domain. Then, the time-frequency plot (representing 10-290 Hz in the frequency domain and 570 ms in the time domain) and amplitude coding plot with size 128x128 are generated.

[0023] FIG. 6B is a plot of an empirical analysis of the model input size, network performance, and model complexity: While keeping the information resolution (2.18 Hz / pix and 4.46 ms / pix) the same, we trained and evaluated the two classifiers with different input sizes from 32x32 (10-80 Hz, ±72 ms) to 224x224 (10-500 Hz, ±500 ms), we plotted the error rate of these two classifiers and the average run time with the corresponding input size together. The error rate of these two classifiers was defined as 1 - 0.5 ( Accartifacts+ AccspkHFo ) in five-fold cross-validation and the average run time (T) was the time to predict J1,000 samples on CPU using a Linux Machine 10 times, — ' Tartifacts+ TspkHFO We chose 128 as the input size because it gave us the best tradeoff between speed and performance. (Acc: accuracy)

[0024] FIG. 7 is a schematic showing that data augmentation was dynamically conducted during neural network training. Within each training epoch, batches of data were sampled from the dataset. For each data sample (HFO event), a time-domain augmentation shift, t, was uniformly sampled from -50 ms to +50 ms. The region of interest for generating features (centered at the middle of the HFO with a 285 ms window) was then shifted on the time axis to select an EEG tracing; for example, in this figure, when t=+50 ms, the window became -235 ms to 335 ms from -285 ms to 285 ms. Subsequently, an integer (0 or 1) was randomly chosen to determine whether the EEG tracing should be flipped in the time domain, with 0 indicating no flipping and 1 indicating flipping. The processed EEG tracing was used to generate two feature images: the time-frequency plot and the amplitude coding plot (the latter only for spk- HFO classification). These images were then fed into the CNN for neural network training

[0025] FIG. 8 is a flowchart of the training and pruning procedure; after the training completes, the pruning and fine-tuning are conducted iteratively to reduce the model size to save computational cost.

[0026] FIG. 9A is a schematic of the spatial distribution of intracranial electrodes. Electrode contact locations within the standardized MNI brain space from various views (Lateral left, Lateral right, Anterior, Posterior).

[0027] FIG. 9B is a schematic of the Variational Autoencoder (VAE) training. Time-frequency EEG data representing HFOs serve as input into the VAE, which outputs a reconstructed image. The VAE's encoder generates a latent distribution of mean and variance while the decoder reconstructs the time-frequency plot from a sampled latent vector. The loss function is a combination of perceptual loss (to capture morphological differences) and KL divergence (measuring the latent distribution's deviation from a normal distribution).

[0028] FIG. 9C is a schematic of the HFO classification pipeline: A two-stage, morphologybased classification process uses Gaussian Mixture Models (GMMs) for unsupervised learning. The first stage identifies artifacts (mArtifact) by latent codes and reconstruction loss; the second stage distinguishes pathological HFOs (mpHFOs) cluster; the cluster with a higher resection percentage in seizure-free patients after resection was deemed pathological. This process trains two unsupervised classifiers to be used on the test set.

[0029] FIG. 9D is a schematic of the interpretability analysis. Left: Latent space clustering, visualized viat-SNE, groups clinically relevant HFO classes, facilitating interpretation. Right: Dimensional perturbations in latent space discover specific neurophysiological features inreconstructed images, providing a bridge between clinical insights and the VAE's feature representations.

[0030] FIG. 10 is a schematic of the model training workflow, according to aspects of the present disclosure.

[0031] FIG. 11 is a schematic of the model interpretation analysis workflow, according to aspects of the present disclosure.

[0032] FIG. 12 is a schematic of the building postoperative seizure outcome prediction model using HFO information.

[0033] FIG. 13 is a serios of plots of an ablation study on different latent dimension. The left panel shows the distribution of dimensions for a VAE model trained with a latent dimension of 8, across five folds. The right panel displays the distribution for the same model configuration but with a latent dimension of 16. Several dimensions in the 16-dimension model do not spread well, indicating that these dimensions are underutilized during training. This compactness suggests a potential reduction in dimensionality. Based on these observations and the more effective utilization of dimensions seen in the boxplot for the 8-dimension model, we chose to proceed with a latent dimension of 8, as it demonstrates optimal use of all available dimensions during training.

[0034] FIG. 14A is an AUCs (area under the curve) of logistic regression (a single-variable classifier) using resection ratios as a variable to predict postoperative seizure freedom are shown based on different types of HFO resection ratios: HFO (unclassified HFO detection), spkHFO (HFO with spike), and mpHFO (morphology-based pathological HFO defined by VAE).

[0035] FIG. 14B is an AUCs from multivariable logistic regression models with different types of HFO resection ratios as the main predictive variable while incorporating subject-wise demographic information (baseline demographics: age and gender) are shown. As an additional predictive variable, the resection status of SOZ was also included. Note that combining the baseline demographics, SOZ resection status, and resection ratio of mpHFO provided the most favorable classification performance (base+soz+mpHFO, AUC=0.71).

[0036] FIG. 14C is a plot of the mean Fl scores (Fl) and standard deviation for the random forest models, trained on the training set using five-fold cross-validation and evaluated on the test set across the five folds and using different features. Note that the resection ratio of mpHFO has a comparable predictive performance compared to the current clinical standard (SOZ resection status). By combining all the features (baseline demographics, mpHFO resectionratio, and SOZ resection status), the model achieved high performance with an Fl score of 0.797 ± 0.081.

[0037] FIG. 15A shows the projected latent space, color-coded by predicted results (mpHFO, non-mpHFO, and mArtifact) from VAE-based HFO morphology inference pipeline on one example fold, shows the 2D projected eight-dimensional latent codes by the t-SNE.

[0038] FIG. 15B shows the same latent space from FIG. 15 A color-coded by conventional HFO classification (spkHFO, non-spkHFO, and artifact).

[0039] FIG. 15C is a box plot of the HFO rates (number of HFO detections / channel / minute) of mpHFO and non-mpHFO based on the location (SOZ versus non-SOZ) across three datasets (outliers were removed for better visualization quality). The rate of mpHFO (mpHFO / channel / min) was higher in SOZ than in non-SOZ. The rates of non-eHFOs (non- eHFOs / min / channel) did not differ between the SOZ and non-SOZ.

[0040] FIG. 15D is a morphological analysis of the time-frequency plot for an example subject. The pathological counterparts (mpHFOs) have higher values throughout the HFO band (> 80 Hz), around the center point (0 ms, where HFOs were detected) than non-mpHFOs; furthermore, higher values of mpHFOs at the non-HFO band (10-80 Hz) throughout the time window compared to non-mpHFOs are exhibited.

[0041] FIG. 15E is a morphological analysis of the time-frequency plot for the overall subjects. The overall template (mean) of all subjects resembles a “hanging bell” shape (mean > 0.5).

[0042] FIG. 15F is a plot of FFT amplitude as a function of frequency from an example patient (shaded areas indicate one SD). Note the FFT peak (around 80-90 Hz) is seen in non-mpHFOs but not in mpHFOs. mpHFO has a higher sub-HFO frequency signal intensity, so the HFO peak may be masked.

[0043] FIG. 15G is a plot of FFT amplitude ratio of (mpHFO - non-mpHFO I non-mpHFO) as a function of frequency, demonstrating the presence of relatively high FFT amplitude in mpHFO within the low frequent band, peaked around 23 Hz.

[0044] FIG. 15H is a time-frequency plot of an example predicted mpHFO.

[0045] FIG. 151 is an EEG tracing of the same mpHFO of FIG. 15H with the detected HFO part colored in red.

[0046] FIG. 15 J is an EEG tracing bandpassed between 10 and 80 Hz of the same mpHFO of FIG. 15H with HFO detection colored in red.

[0047] FIG. 15K is a time-frequency plot of an example predicted non-mpHFO.

[0048] FIG. 15L is an EEG tracing of the same non-mpHFO of FIG. 15K. with the detected HFO part colored in red.

[0049] FIG. 15M is an EEG tracing bandpassed between 10 and 80 Hz of the same non- mpHFO of FIG. 15K with HFO detection colored in red.

[0050] FIG. 16A is a visualization of the latent space for HFOs, color-coded by different datasets (UCLA SEEG, UCLA grid / stnp, and Detroit grid / strip) on a specific fold.

[0051] FIG. 16B is a five-fold cross statistical test showing that the classifiers trained using the actual dataset label (Real) did not show significantly better accuracy (Acc) than those trained by using the permuted label (Random).

[0052] FIG. 16C is a visualization of the latent space for HFOs color-coded by different genders (male and female) on a specific fold.

[0053] FIG. 16D is a five-fold cross statistical test showing that the classifiers trained using the actual genders label (Real) do not show significantly better Accuracy (Acc) than those trained by permuted label (Random).

[0054] FIG. 16E is a visualization of the latent space for HFOs color-coded by different age groups (0-5, 6-10, 11-15, 16-20, and 21+) on a specific fold.

[0055] FIG. 16F is a five-fold cross statistical test showing that the classifiers trained using the actual age group label (Real) do not show significantly better Accuracy (Acc) than those trained by label permuted data (Random).

[0056] FIG. 16G is a visualization of the latent space for HFOs color-coded by different pathologies (HS, FCD, Tumor, and Others) on a specific fold.

[0057] FIG. 16H is a five-fold cross statistical test showing that the classifiers trained using actual pathology labels (Real) do not show significantly better Accuracy (Acc) than those trained by label permuted data (Random).

[0058] FIG. 17 A is a visualization of the latent space for HFOs from preserved brain regions of post-surgical, seizure-free subjects, color-coded by anatomical locations (frontal, temporal, parietal, occipital, and limbic regions) on a specific fold.

[0059] FIG. 17B is a five-fold cross statistical test showing that the classifiers trained using actual anatomical location label (Real) show significantly better accuracy (Acc) than those trained by permuted label (Random) (Acc. Real mean=0.286, std= 0.046, Acc. Random mean = 0.230, std= 0.033. p = 6.67e-10).

[0060] FIG. 17C is an averaged confusion matrix on the test set across five trials and five-fold (n=25 trials) using actual anatomical locations for HFOs from preserved brain regions ofsubjects who achieved postoperative seizure freedom. Note that HF Os from the occipital region are distinguishable.

[0061] FIG. 17D is a set of averaged time-frequency plots for each anatomical location for HFOs from each brain region. Note that HFOs from the occipital region exhibit distinct features on the time-frequency plot.

[0062] FIG. 17E is a visualization of the latent space for HFOs from SOZ channels, color- coded by anatomical locations.

[0063] FIG. 17F a five-fold cross statistical test showing that for HFOs from SOZ channels, the classifiers trained using the anatomical location label (Real) do not show significantly better accuracy (Acc) than those trained by permuted label (Random) (Acc. Real mean=0.241, std=0.043, Acc. Random mean = 0.23, std=0.027, p = 0.090).

[0064] FIG. 17G is an averaged confusion matrix on the test set across five trials and five-fold (n=25) using actual anatomical locations for HFOs from SOZ channels. Note that HFOs from SOZ are indistinguishable from any anatomical origin.

[0065] FIG. 17H is a set of averaged time-frequency plots for each anatomical location for HFOs from SOZ channels. Note that HFOs from SOZ exhibit similar features within timefrequency plots across the anatomical regions.

[0066] FIG. 18 shows the disentanglement of the latent space to establish neurophysiological characteristics of pathological HFOs (mpHFOs). (Panel a) Pathological dimension visualization. This panel illustrates VAE-identified latent space, which enables separation between mpHFOs and non-mpHFOs. (Panel b) A traversal from the lower to the upper percentile of this dimension reveals morphological evolution from non-mpHFO to mpHFO characteristics for both mpHFO base (orange) and non-mpHFO base (blue) in the decoded output, depicted in the image sequence. (Panel c) Power trend within the pathological dimension at the population level. The box plot aggregates the percentage of power within the designated pathological template (the “hanging bell shape”) region of decoded images, showing an ascending trend with higher values, as indicated by the fitted median line. (Panel d) Model probability scores distribution corresponding to each sample, where the mean of the probability marked as red circles, showed the increased confidence of the model as the value as the increase of the value of that dimension. (Panel e) Peak frequency dimensional visualization. A similar visualization of another dimension is shown. This latent space represents the peak frequency of HFOs. (Panel f) The output of the decoder from the traversal of that dimension, displayed in the image sequence, shows a descending trend in peakfrequency from lower to upper percentiles of the value of that dimension. (Panel g) At the population level, the box plots indicate a negative correlation between the peak frequency dimension value and the peak frequency in decoded images, with a trend line fitted from the median of each box. (Panel h) Model probability scores distribution corresponding to each sample, where the mean of the probability marked as red circles, showed the average confidence of perturbed events was around 0.5. (Panel i) Slow- wave dimension visualization. This latent space represents the slow-wave component of HFOs at 10-20 Hz, separated by mpHFO and non-mpHFO prediction. (Panel j) The output of the decoder from the traversal of that dimension, displayed in the image sequence, shows an increased trend in slow wave power from lower to upper percentiles of the value of that dimension. (Panel k) At the population level, the box plot demonstrates a positive correlation between slow-wave dimension values and slow-wave power in decoded images, with a line fit illustrating the median trend. (Panel 1) Model probability scores distribution corresponding to each sample, where the mean of the probability marked as red circles, showed the average confidence of perturbed events was around 0.5.

[0067] FIG. 19 illustrates the 2D projection of the latent space, color-coded according to classifications from a VAE (mpHFO, non-mpHFO, and mArtifact) and from conventional supervised methods (spkHFO, non-spkHFO, and Artifact), evaluated across five folds. It highlights the alignment between VAE predictions and conventional supervised classifications, demonstrating that classifiers developed from unsupervised learning can achieve performance comparable to those developed from supervised methods.

[0068] FIG. 20 is a morphological analysis of mpHFOs' time-frequency plot is shown based on subgroups, each corresponding to a specific variable. Consistent morphological features of mpHFOs are seen across the subgroups. mpHFOs have higher values throughout the HFO band (> 80 Hz), around the center point (0 ms, where HFOs were detected) than non-mpHFOs; furthermore, higher values of mpHFOs at the non-HFO band (10-80 Hz) throughout the time window compared to non-mpHFOs are seen.

[0069] FIG. 21 is a set of plot displaying multiple aspects of HFOs: the HFO power ratio in the first column, HFO peak frequency in the second column, STE HFO length in the third column, and MNI HFO length in the fourth column. The HFO power ratio is defined as the ratio of the signal power within ±50ms of the HFO center to the power within the entire time window of ±285ms. The peak frequency is identified as the frequency component that has the most power within the 80-300 Hz range. Both STE and MNI HFO lengths are measured inmilliseconds as the duration of the HFO event. Each histogram's mean value is annotated within its respective panel. Additionally, statistical comparisons revealed that physiological HFOs exhibited lower power (p<0.001, one-sided t-test), higher peak frequency (p<0.001, one-sided t-test), and shorter durations in both STE and MNI HFOs (p<0.001 for all, except in the frontal location where p<0.01, one-sided t-test) compared to those in other anatomical locations.

[0070] FIG. 22 is a set of plots presenting the distribution of each latent dimension, trained within five folds and depicted by separating mpHFO (in orange) and non-mpHFO (in blue). It is evident that mpHFO and non-mpHFO exhibit significant separation in some dimensions, while others show less distinction.

[0071] FIG. 23 is a sequence of synchronization networks. A one second long iEEG segment is used to construct one synchronization network. Assuming the interictal segment is T seconds long and non-overlapping one second long sliding window approach is used, we have a sequence of T synchronization networks.

[0072] FIG. 24 is a synchronization network workflow. The first step involves selecting a frequency band of interest to filter the iEEG segment and transform it into the frequency domain. Once we have the frequency domain representation, then we construct a synchronization network to capture the interictal interactions between the iEEG channels.

[0073] FIG. 25 is a schematic of frequency bands. Epilepsy research has focused on specific frequency bands depending on the specific research question. Herein, the synchronization in High Gamma, Ripple and Fast Ripple band are studied.

[0074] FIG. 26 is a schematic of edge directionality. If the total phase in channel j leads the total phase in channel i then we assign a directionality from i to j and if the total phase in channel j lags the total phase in channel i then we assign a directionality from j to i.

[0075] FIG. 27 is a schematic of SOZ classifier setup. The SOZ classifier takes in as input a node feature constructed from the interictal network property dynamics and outputs a probability of the node being an SOZ.

[0076] FIG. 28 is a network of a hypersynchronous state with SOZ node P3 as hub. A snapshot of the synchronization network constructed from a one second long interictal segment with blue nodes representing channels that drove seizures in hospital monitoring, red nodes representing channels that did not drive seizures in hospital monitoring but were resected by the doctors, and green nodes representing channels that neither drove seizures in hospital monitoring nor resected by the doctors. The emergence of a hypersynchronous state can be observed, characterized by a star like network, with a SOZ node P3 as the hub.

[0077] FIG. 29 is a network of a hypersynchronous state with SOZ node G53 as hub. Another snapshot of the synchronization network constructed from a different one second long interictal segment with blue nodes representing channels that drove seizures in hospital monitoring, red nodes representing channels that did not drive seizures in hospital monitoring but were resected by the doctors, and green nodes representing channels that neither drove seizures in hospital monitoring nor resected by the doctors. The emergence of a hypersynchronous state can be observed, characterized by a star like network, with a SOZ node G53 as the hub.

[0078] FIG. 30 is a network of a hypersynchronous state with Resected Non-SOZ node Pl as hub. Another snapshot of the synchronization network constructed from a different one second long interictal segment with blue nodes representing channels that drove seizures in hospital monitoring, red nodes representing channels that did not drive seizures in hospital monitoring but were resected by the doctors, and green nodes representing channels that neither drove seizures in hospital monitoring nor resected by the doctors. The emergence of a hypersynchronous state can be observed, characterized by a star like network, with a resected Non-SOZ node Pl as the hub.

[0079] FIG. 31 is a plot of the sharp rise in Fiedler eigenvalue during episodes of large scale synchronization: The Fiedler eigenvalue is graphed over time throughout the interictal period. Examination of the plot reveals distinct spikes in the Fiedler value at certain intervals. These intervals align with periods of hypersynchronous interictal network states, thereby confirming the efficacy of the Fiedler value in capturing such states.

[0080] FIG. 32 is a plot showing the time evolution of out-degree of a node in the network. From the plot, it can be observed that out-degree occasionally enters into an outlier zone marked in red.

[0081] FIG. 33 is a plot of the time evolution of eigencentrahty of a node in the network. From the plot, it can be observed that the eigencentrahty occasionally enters into an outlier zone marked in red.

[0082] FIG. 34 is a flow chart of a 20 dimensional feature vector per node for given frequency band. For a given frequency band, each node is represented by a 5 dimensional feature vector that is evolving over time. The time evolution of each feature is summarized through 4 sufficient statistics thus giving a 20 dimensional feature vector per node per frequency band.

[0083] FIG. 35 is a node feature vector construction workflow: Synchronization effects were studied in three frequency bands: high gamma, ripple and fast ripple. Each frequency band results in a sequence of synchronization networks which in turn gives a 20 dimensional nodefeature. Concatenating the features across the 3 frequency bands we get a 60 dimensional feature vector per node.

[0084] FIG. 36 is a schematic of node labeling: Only seizure free patients were considered for training and for those patients SOZ nodes were labeled as 1 and Preserved Non-SOZ nodes as 0.

[0085] FIG. 37 is a schematic of the training and testing corpus for 17 UCLA seizure free patients. For each fold, there were 16 UCLA seizure free patients in the training corpus and the left out UCLA seizure free patient in the testing corpus resulting in one random forest model per fold.

[0086] FIG. 38 is a schematic of the surgical outcome prediction model. The surgical outcome prediction model takes feature vector of a patient as input and outputs a probability of the patient becoming seizure free which can be then thresholded to predict whether the patient will become seizure free or not.

[0087] FIG. 39 is a plot of the probability of seizure freedom. A box plot illustrating the probability of seizure freedom, as predicted by the 3-dimensional model, stratified by the actual surgical outcomes. Mean seizure freedom probability of seizure free group is slightly higher than the non-seizure free group.

[0088] FIG. 40 is a plot of the 2D visualization of the feature vector for 159 patients. resSOZ ratio is plotted as a function of resection ratio for all 159 patients. The patients are colored based on ground truth surgical outcome: green circles represent seizure free patients and red crosses represent non-seizure free patients. There is a subset of seizure free patients (circled), that has a low resection ratio with some predicted SOZ channels in the preserved region. These patients are incorrectly predicted as non-seizure free by the 3 dimensional model. A majority of the false positives, in the dashed boxed, are situated in the far right comer of the plot with a very high resection ratio and all the predicted SOZ channels removed.

[0089] FIG. 41 is a flow analysis pipeline. The set of nodes in the interictal synchronization network are partitioned into 3 non-overlapping set based on the interictal SOZ classifier predictions and resection information. The flow between these 3 regions is computed by summing the weights of the undirected edges connecting the nodes in these regions.

[0090] FIG. 42 is a plot (left) and schematic (right) of the undirected flow dynamics in a seizure free patient. The flows between the three regions (fl, f2, f3) is plotted as a function of time where each time index corresponds to one interictal synchronization network. On average,there is minimal flow between preserved SOZ and preserved Non-SOZ given the resected SOZ cut.

[0091] FIG. 43 is a plot (left) and schematic (right) of the undirected flow dynamics in a non- seizure free patient. The flows between the three regions (fl, f2, f3) is plotted as a function of time where each time index corresponds to one interictal synchronization network. On average, there is moderately strong flow between preserved SOZ and preserved Non-SOZ given the resected SOZ cut.

[0092] FIG. 44 is a plot of the probability of seizure freedom (5 dimensional). A box plot illustrating the probability of seizure freedom, as predicted by the 5 dimensional model, stratified by the actual surgical outcomes. Mean seizure freedom probability of seizure free group (0.82) is higher than the non-seizure free group (0.48).

[0093] FIG. 45 is a schematic of suggested resection margin based on a workflow for designing minimum resection margins for seizure freedom.

[0094] FIG. 46 is a set of plots of ‘mini seizures’ in interictal period. The Fiedler eigenvalue is plotted over time during the interictal period (left), with the raw waveform amplitude of four channels over the same time period (right). Notably, time segments characterized by hypersynchronous activity (sharp rise in Fiedler value) are boxed in black. Observing the amplitude plot, it becomes evident that the channels exhibit excessive synchronization in the time domain during these hypersynchronous states, bearing resemblance to that observed during seizures.

[0095] FIG. 47 is a set of plot of the correlation between interictal network properties and channel level pathological events. The spk-HFO rate of a channel is graphed over time during the interictal period (left), with the eigenvector centrality of the same channel over the same time period (right). Upon examination of both plots, it becomes apparent that the spk-HFO rate of a channel appears to be correlated with eigenvector centrality of that channel.

[0096] FIG. 48 is a plot of the probability of seizure freedom (7 dimensional): A box plot illustrating the probability of seizure freedom, as predicted by the 7 dimensional model, stratified by the actual surgical outcomes. Mean seizure freedom probability of seizure free group (0.9) is significantly higher than the non-seizure free group (0.27).

[0097] FIG. 49 is a plot of the Pattern in the correlation feature space: The ratio of variance to7ar[Rc] mean is plotted against the mean Ff^c] for all 159 patients. Patients are color-coded £ty2] based on their ground truth surgical outcome: green circles depict seizure free patients, whilered crosses represent non-seizure free patients. A distinct separation between the two patient groups is evident in this two-dimensional feature space.

[0098] FIG. 50 is a flowchart of the k-nearest neighbor clustering in the correlation feature space.

[0099] FIG. 51 is a plot of four clusters in the correlation feature space. In the 2-dimensional correlation feature space, four distinct clusters emerge. Among them, Cluster 2 and Cluster 4 are primarily composed of seizure free patients, whereas Cluster 1 and Cluster 3 predominantly consist of non-seizure free patients.

[0100] FIG. 52 is a box plot of channel level R2values of a patient from Cluster 1 : For a patient in cluster 1, channel level R2values categorized by channel resection status are plotted, with green points indicating resected channels and red points representing preserved channels. A significant number of channels with spk-HFO detections, have a low R2value, suggesting that a sizeable number of detected spk-HFOs remain unexplained by the sampled interictal network.

[0101] FIG. 53 is a box plot of channel level R2values of a patient from Cluster 4: For a patient in cluster 4, channel level R2values categorized by channel resection status are plotted, with green points indicating resected channels and red points representing preserved channels. Majority of the channels with spk-HFO detections, have a high J?2value, suggesting that most of the spk-HFOs detected are explained by the sampled interictal network.

[0102] FIG. 54 is a box plot of channel level R2values of a patient from Cluster 2: For a patient in cluster 2, channel level R2values categorized by channel resection status are plotted, with green points indicating resected channels and red points representing preserved channels. The majority of channels exhibit zero spk-HFO detections, resulting in zero R2values. Among channels with a few spk-HFO detections, their R2values are high, suggesting that the spk-HFO rates of these channels are explained by the sampled interictal network property.

[0103] FIG. 55 is a box plot of channel level R2values of a patient from Cluster 3 : For a patient in cluster 3, channel level R2values categorized by channel resection status are plotted, with green points indicating resected channels and red points representing preserved channels. The majority of channels exhibit zero spk-HFO detections, resulting in zero R2values. Among channels with a few spk-HFO detections, their R2values are low, suggesting that the spk-HFO rates of these channels remain unexplained by the sampled interictal network property.Detailed Description

[0104] The systems and methods described herein represent a framework for analyzing brain data including a plurality of modules. FIG. 1 illustrates a non-limiting example of theframework 100 is a system comprising multiple modules and sources of data. The framework 100 includes one or more multimodal integrated clinical datasets 102. As described further in the description of FIG. IB, the multimodal integrated clinical dataset may be a database that stores multimodal invasive and non-invasive brain scan data, including but not limited to magnetic resonance imaging (MRI), functional MRI (fMRI), transcranial magnetic stimulation (TMS), positron emission tomography (PET), scalp electroencephalography (EEG), and intracranial EEG (iEEG). In anon-limiting epilepsy example, data are collected during seizure- free or interictal periods. Alternatively, the data are collected during seizure onset and / or ictal periods. Observation periods are usually limited to visits to a medical facility and seizures are not guaranteed to occur during any particular visit; moreover, even if seizures do occur it might be only due to only a few of the seizure onset zones and other such regions might lie inactive only to trigger seizures at a later time. Epileptogenic tissues are active even when there is no observable seizure, and such interictal activities reveal biomarkers in signals measured using different modalities. The multimodal clinical dataset 102 may also include patient metadata 103 and physician diagnostic inputs. In a non-limiting example, this dataset 102 serves as the central repository for all data used in the framework 100.

[0105] The framework 100 further includes an Al model generator 104. In a non-limiting example, the Al model generator 104 receives data from the one or more multimodal integrated clinical datasets 102 and generate Al models 105 to detect and classify biomarkers. In a nonlimiting example, this may involve implementing machine learning (ML) algorithms, deep learning (DL) models, or other advanced analytics techniques. FIG. IB also provided further details related to the generating the Al models 105. In a non-limiting example, the biomarkers include high frequency oscillations (HFOs) from scalp EEG or iEEG. Alternatively, the biomarkers may be spike, sleep spindle, interictal epileptiform discharges (lEDs), K- complexes, sharp wave-ripples (SWRs), cross-frequency coupling (CFC), and slow oscillations. In another example, the biomarker may be a network-based biomarker (as described in further detail in Example 3 below).

[0106] The generated Al models 105 are further trained in an Al training module 106. Specifically, the Al training module 106 peforms a DL-based classification of one or more biomarkers based on one or more features extracted from the biomarkers. In a non-limiting example, the biomarkers may be classified via a supervised DL classification, self-supervised DL classification, or an unsupervised DL classification. In another non-limiting example, this step could use both supervised and unsupervised machine learning (ML) methodologies to trainmodels that would discover and determine multimodal signal patterns or other features for both epileptogenic and seizure onset zones, and locations of various cognitive tasks in the brain. In a supervised setting, patterns that have been conjectured by doctors to be relevant would act as labels for model training. In a weakly supervised scheme, patient outcome data can be utilized to further augment labels; for example, signals from regions when resected lead to postoperative seizure freedom can be considered as carriers of epileptogenic biomarkers. In a completely unsupervised framework, generative models such as variational autoencoders (VAEs) can be used and the different clusters in the latent space can be used as putative biomarkers patterns in the underlying scans. In a non-limiting example for training the Al models 105, the framework 100 uses (i) partially labeled data annotated by clinicians, (ii) creates weakly labeled data based on treatment outcome, and (iii) self-supervised generative models to create putative labels which are then verified based on treatment outcomes and inputs from experts.

[0107] In a non-limiting example, data input and trained in the Al model generator 104 and Al training module 106 may be annotated and verified by a medical expert and incorporated throughout the model development process. Further details regarding the training process are provided in the examples below. The performance of the Al models 105 may be assessed by comparing their predictions to known outcomes, such as post-operative seizure freedom or physician-identified biomarkers. The Al models 105 may be refined based on this feedback, incorporating additional user input and surgical outcomes as necessary. For example, accuracy of such models in predicting eloquent cortex can be verified by providing stimulations to the predicted regions and logging response data. False positive and false negative predictions can then be used as additional labels to refine the models. Similarly, predictions for epileptogenic zones can be correlated with patient treatment outcomes leading to more data about false positives and negatives.

[0108] The framework 100 further includes input from multimodal individual clinical patient data 108. In a non-limiting example, the input may be input data of a patient’s brain comprising one or more modalities to identify the one or more features in the input data. The clinical patient data 108 may include data of any of the modalities previously listed for the multimodal integrated clinical data 102. In a non-limiting example, clinical patient data may also include patient metadata 109. A patient’s recorded data is input to a model including the model(s) previously generated and trained by the Al Model Generator 104 and Al training module 106 to identify biomarkers and their associated features therein.

[0109] Upon identification of these biomarkers, this data and the trained Al models are aggregated into an Al-based clinical decision recommendation module 110. The module 110 may generate a map, such as a likelihood heatmap of the patient’s brain, indicating epileptogenic zones and functional activity (eloquent regions). Additionally, or alternatively, the likelihood heatmap may visualize connectivity, metabolic activity, neurotransmitter systems, receptor distribution, gene expression, and disease-related changes. The module 110 may further generate a treatment plan for one or more target regions in the patient’ s brain based on the generated brain map. A non-limiting example of the data aggregation into the recommendation module 110 is shown in FIG. 1C.

[0110] The framework 100 further includes a user interface 112 configured to provide the treatment plan, information related to the one or more biomarkers, cognitive spatial biomarkers, and other statistics and insights to a user for evaluation 114. This output may be communicated to a user visually on a display screen. The user interface 112 may further be configured to receive input by a user, such as the annotated data or other parameter settings such as annotation window length or processing parameters in filtering etc. This input may be communicated via computer components such as a mouse, keyboard, joystick, etc.[OHl] In a non-limiting example, the framework further includes performing the recommended treatment 116 from the treatment plan and obtaining patient treatment response data 118. The treatment may involve electrical stimulation to the target regions or surgical recission of the target regions. Clinical outcome data may then be measured for the patient post-treatment and provided back into the framework for refining the Al models. For example, outcome data may include post-operative seizure cessation or user-identified biomarkers. In one example, the outcome data may be incorporated into the multimodal integrated clinical dataset 102. Alternatively, the outcome data may be entered by a used via the user interface 112. For example, the outcome data may refine at least one of the one or more biomarkers, the cognitive spatial biomarkers, the brain map, or the treatment plan. This process represents a hypothesis-driven iterative approach where an expert-assisted Al tool makes recommendations based on the existing knowledge and data that are adopted for treatment and surgical interventions, which in turn generates post-operative and post-treatment outcome data that is then used to adapt and modify recommended practices. In particular, predictions made by the inference engine are compared with patient outcomes to determine false positives and false negatives so that a Bayesian approach can be used to update the models as described in the above. By continuously updating the models with new patient data and outcomes, theframework will be better equipped to provide more precise and effective recommendations for future patients.

[0112] Details of generating the treatment plan combining pre-trained Al models 105’ and clinical patient data 108 are shown in FIG. 1C. When a patient's data 108 is recorded from multiple sources, the data is sent to the system, which then uses the corresponding metadata 109 and pre-trained models 105’ to evaluate the patient's data for each specific source. The pretrained models 105’, the patient data 108 and patient metadata 109 are aggregated in an aggregation module 111 and integrated into the Al-based clinical decision recommendation module 110, which and makes a final clinical decision recommendation or provides a treatment plan 113. In one embodiment, the Al module’s output is a 3D brain mapping 115 with likelihood heatmaps that visualize complex data, such as functional activation (eloquent regions), connectivity, metabolic activity, neurotransmitter systems, receptor distribution, gene expression, and disease-related changes such as epileptogenic zone.

[0113] Referring to FIG. 2A, a method 200 for identifying pathological biomarkers is shown. The method 200 includes generating the Al model, training the Al model, and refining the Al model. In a non-limiting example, the method 200 may be carried out using a system as previously described by the framework 100. Step 202 includes receiving, by the system comprising a processor, one or more datasets from one or more modalities associated with brain activity of a population. In a non-limiting example, the datasets and modalities may those provided above. Next, the processor identifies one or more biomarkers for each of the one or more datasets associated with brain activity of a population at step 204. At step 206, the processor may extract one or more features from the one or more biomarkers. The process may also identify cognitive spatial biomarkers for each of the one or more datasets at step 208. The processor then performs DL-based classification of the one or more biomarkers based on the one or more features at step 210. At step 212, a user interface may output at least one of the one or more biomarkers or cognitive spatial biomarkers identified in the data. Optionally, at step 214, the processor may receive outcome data based on the one or more biomarkers or cognitive spatial biomarkers. As described previously, predictions of the eloquent cortex and epileptogenic zones may be verified by providing stimulations to the predicted regions. The resulting outcome data of the stimulations may be used to refine the Al model generation or training. For example, the outcome data may refine the biomarker identification at step 204, the feature extraction at step 206, the cognitive spatial biomarker identification at step 208, or any combination thereof.

[0114] FIG. 2B shows a method 220 for identifying treatment target regions of a patient based on pathological features acquired in patient data. A process or defined previously by framework 100 above may be used to perform method 220. At step 222, a processor may receive one or more features extracted from one or more biomarkers associated with seizure activity identified from one or more datasets from one or more modalities associated with brain activity of a population. The processor may also obtain input data of a patient’s brain comprising one or more modalities at step 224. The input data may of any one or more of the modalities described previously. At step 226, the processor identifies the one or more features in the input data. The processor then generates a brain map of the patient based on the one more features identified in the input data at step 228. At step 230, the processor generates a treatment plan for one or more target regions in the patient’s brain based on the brain map. A user interface may output at least one of the one or more biomarkers, the brain map, or the treatment plan.

[0115] In a non-limiting example, the outcome of treatments performed based on the treatment plan generated in 230 may be used as outcome data to further refine the Al model generation and training as indicated by the treatment response data 118 in the framework 100 of FIG. 1A.

[0116] The following examples illustrate examples implementations and embodiments, are not intended to be limiting.Example 1

[0117] A demonstration of the framework 100 may involve refining the HFO biomarker in intracranial EEG (FIGS. 3A-3B). The recording device captures intracranial EEG data from the patient and uploads it to the cloud, along with basic patient information. In this example, the framework is implemented using Python, PyQt, and Py Torch. The framework includes cloud-sdie and client-side services having the following capabilities:

[0118] 1. The client-side software application is responsible for receiving and visualizing the EEG data. To prepare the data for further analysis, the application performs several processing tasks, such as filtering, time-frequency analysis, bipolar montage visualization, and resampling.

[0119] 2. After visualizing and processing the EEG data, the client-side application proceeds to detect potential HFO events using re-implemented HFO detection algorithms in Python, including the Short Term Energy (STE) and Montreal Neurological Institute and Hospital (MNI) algorithm. To optimize the computational pipeline, the signal-filtering and detection algorithms are parallelized in each channel to each CPU core, utilizing a multi-processing mechanism. This optimization significantly reduces the runtime of event detection. The resultis a more efficient and accurate HFO event detection process. The HFO detection pipeline is depicted in FIG. 3A.

[0120] 3. Once HFO events have been detected, the client-side application extracts features from these events. The feature extraction procedure involves generating two images for each event: a time-frequency plot (using the Morlet wavelet transform) and an amplitude coding plot (a 2D representation of the signal's intensity) 500 ms around the middle point of each HFO event. The feature extraction procedure is heavily optimized by adopting multi-processing for each event and replacing the iteration (for loop) with matrix multiplication when computing the time-frequency plot. This optimization reduces the time required to extract features from each HFO event and increases the overall efficiency of the analysis.

[0121] 4. After features are extracted from each HFO event, the client-side application utilizes two pre-trained deep-learning models, implemented with PyTorch, that have been trained on the cloud-side service. These models classify each HFO event by feeding the corresponding features into the model. In this example, the two models provided are designed for the most demanding HFO refinement tasks needed for clinical decision-making: artifact rejection and spike-HFO classification. The detected results are displayed in the graphical user interface (GUI) of the client-side software application, as shown in FIG. 4A-4B of the software platform.

[0122] 5. If required, the doctor can annotate additional HFO events and upload them to the cloud-based service for further neural network training and optimization. The annotation procedure is done by displaying the original signal, filtered signal, and the time-frequency plot of the 500 ms signal around the center of the event, as shown in FIGS. 4A-4B of the software platform. This annotation process provides additional data that can be used to train and optimize the deep-leaming models used in the analysis, resulting in improved accuracy and efficiency in future analyses.

[0123] The cloud service can be hosted on any web service, offering several functionalities for the software platform.

[0124] 1. Deep Learning and Machine Learning Infrastructure: Cloud side provides the necessary infrastructure for deep learning and machine learning, including PyTorch as the deep learning framework and GPU resources for supporting the training of neural networks.

[0125] 2. Training, Pruning, and Fine-tuning: The cloud-side service carries out the training, pruning, and fine-tuning of the neural network, ensuring an efficient model is deployed back to the client-side application. In this example, the network architecture is ResNet-18, with 15 training epochs, using 5-fold cross-validation on a private dataset containing 19 patients' 10-minute EEG data and 12,494 events with corresponding expert annotations. The model pruning is conducted using DepGraph, with 5,000 pruning iterations, and the model is trained for 5 epochs at every 500 pruning iterations. The task is a classification task, so the loss function is binary cross-entropy loss. The performance metrics used to evaluate the network are accuracy, Fl-score, precision, and recall.

[0126] 3. Data Storage: The cloud-based service stores the integrated dataset containing the EEG data, as well as other clinical data such as fMRI results and additional patient-specific information. This storage capability ensures that all relevant data is easily accessible and securely managed within the cloud.Example 2Method

[0127] PyHFO is a multi-window Graphical User Interface (GUI) desktop application specifically designed for the efficient analysis and classification of HFOs. It presents a user- friendly and intuitive interface that caters to both technical and non-technical users, streamlining the process of HFO detection and classification. PyHFO operates through four primary stages upon loading an EEG recording: EEG signal reading, data filtering, HFO detection, and Deep Learning (DL)-based HFO classification. These stages are detailed in a data flowchart, as seen in Figure 1. The output of this pipeline includes detected events based on the implemented detection algorithm, accompanied by annotations of real HFOs, artifacts, and spkHFO, generated using pre-trained DL-based HFO classification models. The specifics of each critical stage are elaborated upon in the subsequent sections.

[0128] HFO detection algorithms

[0129] In PyHFO, we implemented two automatic HFO detection algorithms into the standalone executable software. We selected the Short Time Energy (STE) and The Montreal Neurological Institute and Hospital detector (MNI) to implement because they are the two primary detection algorithms in the widely used Matlab-based HFO analysis tool, RIPPLELAB. where they have demonstrated success in numerous studies. The PyHFO’s modular architecture, coupled with the open-source property of the project, allows for easy integration of alternative HFO detection methods if equired. Developers need to adhere to a straightforward interface. Moreover, any added methods will inherently benefit from the established multi-processing paradigm.

[0130] We have faithfully replicated the precise parameters and computational implementation of both algorithms from RIPPLELAB in Python presented in Table 2.1 and Table 2.2.However, we have introduced certain modifications. This includes replacing functions, such as the gamma distribution parameter estimation, with the official Scipy’s APIs. Additionally, we have exposed the random seed to the user to ensure the reproducibility of the MNI detector. More importantly, to enhance execution efficiency, we have replaced the for loop with matrix multiplication, particularly in the Gabor wavelet computation. These alterations may result in slightly divergent detection outcomes compared to those obtained from RIPPLELAB. A comprehensive analysis and comparison of these results will be presented in the dedicated analysis section.Table 2.1. Detailed parameters with units in parentheses of the STE detector used in UCLA, Zurich, and Rodent datasets: filter_freq: frequency band in bandpass filter; rmswindow: RMS window time; min window: minimum time for an HFO; min_gap: minimum distance time between two HFO candidates; min_osc: minimum oscillations per interval; rms_thres: threshold for RMS in standard deviation (SD); epochjen: duration of determining energy threshold.Table 2.2. Detailed parameters with units in parentheses of the MNI detector used in UCLA, Zurich, and Rodent datasets. Filter freq: frequency band in bandpass filter; epo_CHF: continuous high-frequency epoch; per_CHF: continuous high-frequency percentile threshold; min_win: minimum HFO time; min_gap: minimum distance time between two HFO candidates; thrd_perc: threshold percentile; base_seg: baseline window; base_shift: baseline shift window (representing the ratio of the baseline window is shifted); base_thrd: baseline threshold; base_min: baseline minimum time.UCLA Zurich Rodent filter freq (Hz) [80, 500] [80, 500] [80, 500] epo.CHF (s) 60 60 60 per_CHF (%) 95% 95% 95% min_win (s) 10e-3 10e-3 10e-3 min_gap (s) 10e-3 10e-3 10e-3 thrd perc (%) 99.9999% 99.9999% 99.9999% base.seg (s) 125e-3 125e-3 125e-3 base_shift (0-1) 0.5 0.5 0.5 base thrd (0-1) 0.67 0.67 0.67 base.min (%) 5 5 5

[0131] HFO detector implementation details

[0132] Data reading. PyHFO is designed to accept mainstream EEG data file formats such as the European Data Format (EDF). Additionally, it can process data in the widely-used NumPy format when users employ the deployed Python package (see Software Overview Section). In processing an EDF file, raw data is stored in binary format. Upon reading, to convert the digital (raw) values Draw to the real-world physical voltage V, the digital values are calibrated using maximum and minimum physical voltage V max, V min and maximum and minimum digital values Dmax , Dmin. The equation used by most EEG data processing tools, such as MNE, is given in Equation (1). In this equation, R is the calibration ratio, defined as the ratio of the difference between the maximum and minimum physical values to the difference between the maximum and minimum digital values, and an offset 0 is defined as the difference between the minimum physical value and the product of the calibration ratio and the minimum digital value.

[0133] If s worth noting that RIPPLELAB processes EDF files differently than other mainstream EDF reading tools. Contrastingly, RIPPLELAB performs calibration only through V = R * Draw, with no offset adjustment, resulting in data readings with a DC offset between RIPPLELAB and other EDF reading tools. In our implementation within PyHFO, we have elected to use the more robust calibration equation, Equation (1), as the open-sourced Python package MNE.

[0134] Signal filtering. The voltage value read from

[0135] EDF was then passed through a bandpass filter to extract the signal in the desired frequency domain with the specified ripple and attenuation. The bandpass filter used in the RIPPLELAB was the Chebyshev type II filter; the parameter of this filter consisted of PassBand (FP), StopBand (Fs), PassBand Ripple (rP), and StopBand Attenuation (rs). For constructing such a filter, the order of the filter was first estimated, and then the frequencyresponse was constructed. We noticed that Matlab sometimes could not achieve an exact match to the desired PassBand Ripple and StopBand Attenuation. Therefore, in implementing the PyHFO, we chose to use the filter construction by Scipy as it could produce a more aligned frequency response to the specification. However, since the Scipy could not replicate the filter parameter specified in RIPPLELAB, Fp= 80 Hz, Fs= 500 Hz, rP= 0.5 dB, and rs= 100 dB due to the numerical overflow, we used the closest number in StopBand Attenuation rs= 93 dB instead in our implementation. We visualized such prementioned phenomena in FIG. 5.

[0136] Multi-processing-based detection framework.

[0137] To improve the efficiency of our HFO detection pipeline, we leveraged the multiprocessing capability of Python. In FIGS. 3A-3B, we illustrated how we parallelized the timeconsuming steps, namely data filtering and HFO detection, across each channel. Since the computation for each channel was independent, we assigned each channel’s data to different CPU cores to run simultaneously. This approach maximized the use of available CPU cores, leading to a significant reduction in the detection time. To demonstrate the acceleration of our detector’s running speed, we compared the detection times for MNI and STE detectors using our detector and the RIPPLELAB detector on different hardware machines.

[0138] Lightweight deep learning-based HFO

[0139] classification

[0140] We followed the same artifact and spk-HFO classifier design in Zhang et al. (Refining epileptogenic high-frequency oscillations using deep learning: a reverse engineering approach, Brain Communications, Volume 4, Issue 1, 2022, fcab267) as it had already shown promising performance against expert annotation. The training data was from HFO detected by STE detector via RIPPELAB in UCLA dataset along with annotation from experts (NH and SH). However, several limitations prevent them from being directly used in the natural setting: 1) Currently, the HFO analysis is majorly conducted in CPU machines; access to GPU is not very popular in this domain of study. Application of the proposed network in Zhang in CPU is time-consuming. 2) Then, the generalization ability of models in Zhang cannot be ensured in HFOs detected by other detectors such as MNI. To address 1), we reduced the computational cost of the model by first seeking the smallest information (input size) that can maintain the classification performance and then employing the state-of-the-art neural network pruning technique to reduce the network size. To address limitation 2), we developed a data-augmentation strategy in the neural network training to improve the generalization ability of the model.

[0141] DL model training with data augmentation.Two DL models were trained in PyHFO: the artifact rejection model and the spkHFO classification model. The artifact rejection model classified all events detected from the HFO detector into artifact and real HFO events (the union of spkHFO and non-spkHFO). Meanwhile, the spkHFO classification model classifies the real HFO events into spkHFO and non-spkHFO. The models were evaluated through five-fold cross-validation. The dataset was randomly shuffled and divided into five groups; each group represented the test set of each fold in cross- validation. Within each fold, the remaining 80% of the data was then split into a training set (70% of the whole dataset) and a validation set (10% of the whole dataset). During the training, we used time-domain augmentation to improve the generalization ability of the artifact rejection model and spkHFO classification model. For each event within the training batch, we randomly flipped the EEG signals and randomly shifted the center of the HFO event forw ard and backward 50 ms, as shown in FIGS. 6A-6B and FIG. 7 Both neural networks were trained for 30 epochs; thus, each data sample was augmented 30 times. In the validation and test sets, augmentation was not applied to ensure that events were represented accurately during evaluation. Both networks were trained with a batch size of 128 using Adam optimizer with a learning rate of 0.0003. The final DL models were selected based on achieving the lowest validation loss across epochs during the training process.

[0142] Reduction of the computational cost. The computational cost of a neural network could be measured by the total number of Multiply Accumulate (MACs) for a fixed number of inputs, which was influenced by the input dimension and the architecture size. To reduce the computational complexity introduced in the input dimension, we first reduced the redundancy in the dimension of the input by only using one time-frequency plot for the artifact detector and concatenation of only the time-frequency plot and amplitude coding plot as input for the spk-HFO classifier. Then, we reduced the input dimension from 224x224 to 128x128 by only taking 10 to 290 Hz in the frequency domain and ±285 ms of the center of the event in the time domain; these values were chosen by empirical analysis of balancing the computational complexity of the neural network and the classification accuracy (FIG. 6A-6B). Then we simplified the architecture of the artifact and spk-HFO detector, respectively, by pruning the neural network using DepGraph as shown in FIG. 8. The interactive pruning was conducted for 5,000 iterations, and the model was fine-tuned every 250 iterations by five epochs. Finally, we imposed a rule-based filter by treating all HFOs detected in the beginning one second and last one second as artifacts because the beginning and the ending of the recording would leadto artifact production. The simplified network should run at the best tradeoff between speed and performance in CPU, and we also enabled the use of GPU for users with GPU access on their machines. We evaluated the model complexity by computing MACs using one data sample. Additionally, we measured the average running speed of the model inference on 1,000 data samples to get a more straightforward overview of the model complexity.

[0143] Framework evaluation

[0144] Evaluation patient cohort and iEEG recording. We evaluated the performance of the HFO detectors and classification by using three iEEG datasets.

[0145] UCLA iEEG Dataset (UCLA):

[0146] Intracranial EEG (iEEG) data was obtained via grid / strip electrodes using Nihon Kohden Systems (Neurofax 1100 A, Irvine, California, USA). The recording was acquired with a digital sampling frequency of 2,000 Hz. It contained 19 drug-resistant focal epilepsy subjects. For each subject, separate 10-minute EEG segments from slow-wave sleep were selected at least two hours before or after seizures, before anti-seizure medication tapering, and before cortical stimulation mapping, which typically occurred two days after the implant. This dataset contained 19 ten-minute EEG recording segments across 19 patients with 1,709 monopolar channels (a median of 94 monopolar channels in each recording). The annotation (Artifact, HFO-with-spike, HFO-without- spike) obtained from expert labeling from a previous study (Zhang) of each STE HFO event was also included in this data (see HFO event annotation section for detailed annotation statistics).

[0147] Zurich iEEG HFO Dataset (Zurich):

[0148] Intracranial EEG (iEEG) data (both grid / strip and depth electrode) was obtained with a 4,000 Hz digital sampling frequency with an ATLAS recording system (0.5-1000 Hz passband, Neuralynx) and downsampled to 2,000 Hz. It contained 20 drug-resistant focal epilepsy subjects. Several runs of five-minute EEG segments of interictal slow- wave sleep were recorded for each subject. We followed the same preprocessing procedure in Fedele et al. (Resection of high frequency oscillations predicts seizure outcome in the individual patient. Sci Rep 7, 13836 (2017)) to create bipolar EEG recordings. This dataset contained 385 five- minute EEG recording segments across 20 patients with 9,360 bipolar channels (a median of 23 bipolar channels in each recording).

[0149] UCLA Rodent Dataset (Rodent):

[0150] Rodent intracranial EEG (iEEG) data was obtained with a 2,000 Hz digital sampling frequency with RHD2000 electrophysiology amplifier chips (pass-band between 0. 1 Hz and1,000 Hz). It contained two ten-minute EEG segments from two rodent subjects from the neocortex and hippocampus sampling via depth electrodes. One EEG recording was from a subject with traumatic brain injury (nine channels), and the other one was from a sham-injured control subject (ten channels). Please see the published data in Santana-Gomez et al. (Harmonization of pipeline for detection of HFOs in a rat model of post-traumatic epilepsy in preclinical multicenter study on post-traumatic epileptogenesis. Epilepsy Res. 2019 Oct; 156:106110) for a detailed description of this dataset.

[0151] HFO detector parameters. We used the same parameter settings for RIPPLELAB and PyHFO to compare the consistency of their detection results and runtime because PyHFO essentially replicates the detection pipeline of RIPPLELAB. The STE and MNI detectors utilized identical default parameter settings in RIPPLELAB and PyHFO for the UCLA and Zurich datasets. For the Rodent dataset, we applied the suggested parameters for the STE detector as introduced in the original paper by Santana-Gomez. However, since no suggested parameters existed in the MNI detector in the original paper, we used the default parameter introduced in RIPPLELAB. Tables 2.1 and 2.2 provides exact parameters used in each dataset.

[0152] HFO detector evaluation. To conduct a mathematical evaluation of the detection results between PyHFO and RIPPLELAB, we established a defined representation of events detected by each algorithm. For PyHFO, an event was denoted as (startp , endp ), indicating the exact time location within the EEG recording. Similarly, for RIPPLELAB, an event was represented as (startr, endr). To quantify the degree of overlap between these two sets of events, we introduced the concept of an overlapping ratio, which was defined as

[0153] The resulting value ranged from 0 to 1, with 1 indicating an exact match. To ensure a fair comparison and avoid double counting, we enforced the condition that an event detected by PyHFO can only be matched with a unique event detected by RIPPLELAB. Additionally, the comparison was performed on a channel-by-channel basis, and there was no overlap within events detected by the same detector by the definition of the detecting algorithm. The match of a specific event was defined when the overlapping ratio exceeded a certain threshold, such as 50%. Furthermore, the discrepancy between the two algorithms could be quantified by calculating the ratio of the number of matches to the total number of events detected by RIPPLELAB.

[0154] We conducted four experiments to evaluate the success of our detector implementation, assessing the impact of each module independently in the pipeline. For simplicity, we denoted the data reading (Read), filter design (Filter), and detection algorithm (Algo) of RIPPLELAB. as Readr Filterrand Algor , where the subscript r r represents the RIPPLELAB implementation, and we used the subscript p p to denote the PyHFO implementation. To verify the correctness of our Python implementation, we first extracted the filtered EEG signal from RIPPLELAB. We then fed it into our detectors. (Readr + Filterr + AlgoP), comparing the detection overlap with RIPPLELAB (Readr + Filterr + Algor ), which we referred to as Expl. Since we replicated the logic of the two detectors, i.e., Algor = AlgoP, we expected almost 100% matching between them. To assess the impact of the data reading, we conducted Exp2, replicating the frequency and phase response from RIPPLELAB in Python (Readp + Filterr + Algop ). For Exp3, we evaluated the effect of the filter design by feeding the EEG signal read by RIPPLELAB into our Py thon pipeline (Readr + FilterP+ AlgoP). Finally, we compared the complete implementation of PyHFO (Readp + FilterP+ AlgoP) with the RIPPLELAB implementation (Readr + Filterr + Algor ). We conducted all experiments on our evaluation cohort and evaluated the implemented STE and MNI, respectively. For each experiment, we compared the detected HFOs with those detected from RIPPLELAB on the total number of HFOs, number of exact match HFOs, and the number of at least 50% overlap, evaluating the discrepancies between RIPPLELAB and our implementation in each step of the data processing. By evaluating the effect of each module independently, we were able to demonstrate the success of our detector implementation, providing a comprehensive assessment of the performance of our PyHFO implementation.

[0155] DL-based neural network evaluation. The performance of our trained artifact and spk- HFO detectors was evaluated by comparing the results with the label. We adopted standard metrics for machine learning classification tasks to assess the model’s performance, includingPrecision = - TP+FP , Recall = - TP+FN , Accuracy (Acc) ' = -Tp+TN +pp+FN, and Fl -score =2xPrecision * RecallprepresenlsTmePositives, FP represents False Positives, FN Precision+Recall represents, and TN represents. Given that the model was trained using five-fold cross- validation on STE HFOs, the reported metric values were the mean results of the fivefold cross-validation (on different test sets) with a 95% confidence interval. To evaluate the model performance on the MNI HFOs, experts annotated MNI HFOs from representative patients. Models trained in five-fold cross-validation were used to predict all events. Thereported metric was thus the mean of the metrics from five models, again with a 95% confidence interval.

[0156] Runtime analysis. We conducted our runtime analysis on a Linux machine, a macOS machine, and a Windows machine. The Linux machine had an AMD Ryzen Threadripper 2950X 32-core processor; the Windows machine had an Intel i9-13900K 24-core processor; and the macOS machine had an Apple Ml Pro 8-core processor. For timing the HFO detectors in RIPPLELAB, we followed the modification of RIPPLELAB in Zhang to run the Matlabbased detector and report the runtime for detecting STE and MNI HFOs in the Linux machine. For PyHFO, we ran the detector with the same parameters as in RIPPLELAB and reported the runtime by using single-core (n-jobs = 1) and multi-core (n-jobs = 32 for Linux, n-jobs = 8 for Windows and Mac machines). See Tables 2.1 and 2.2 for exact parameter settings in both detectors. For benchmarking the DL models, we used DL models to predict 1000 samples ten times to get the mean and 95% confidence interval of the runtime in different machines with PyTorch default settings, and we also reported the inference time on an Nvidia Titan RTX GPU for reference.

[0157] Software Overview

[0158] The software version PyHFO is a multi-window GUI developed in PyQt. It is intended to be a user-fnendly and intuitive tool that users with technical and non-technical backgrounds can use to detect and classify HFOs time-efficiently. PyHFO has been released under Academic Licenses (Licenses for Sharing Software Code Non-commercially, UCLA TDG). The GUI interface was implemented in PyQt version 5.15 to make it compatible with hardware platforms such as macOS, Linux, and Windows. The HFO detectors were implemented in Python 3.9.0, and the DL-based detector was implemented in PyTorch 1.6. We chose Python as the programming platform because it is widely used in large-scale data analysis and deep learning in the medical image field. The completed procedure for detecting and analyzing HFOs through PyHFO consists of several steps briefly discussed in FIGS. 3A-3B. We here presented the GUI of the software in FIGS. 4A-4B. After setting the parameters of different detectors and the DL- based classifier, the original EEG signal will be displayed in the main visualization window, and HFOs with different attributes (artifacts, HFO-with-spike, and HFO-without-spike) will be annotated with different colors. The statistics of different kinds of HFOs will be displayed on the summary panel for each channel. All of these statistics can be exported in Excel format for further study by the user.

[0159] Results

[0160] HFO detector evaluation

[0161] In the UCLA dataset, RIPPLELAB detected 12,494 STE HFOs and 10,392 MNI HFOs, while PyHFO detected 12,501 STE HFOs and 10,355 MNI HFOs. In the Zurich dataset, RIPPLELAB detected 31,744 STE HFOs and 70,988 MNI HFOs, while PyHFO detected 31,869 STE HFOs and 70,538 MNI HFOs. In the rodent dataset, RIPPLELAB detected 375 STE HFOs and 42 MNI HFOs, while PyHFO detected 378 STE HFOs and 39 MNI HFOs. In Table 2.3, we demonstrated the breakdown performance of each experiment. Specifically, in Expl, we demonstrated that PyHFO successfully replicated the detection algorithms implemented in RIPPLELAB. The discrepancy in the MNI detector was due to different random seed mechanisms RIPPLELAB and PyHFO used. Controlled variable experiments showed that different data readings (Exp2) and filter (Exp3) do affect the performance of the detector but with a minimal effect of around 3 to 7% difference between RIPPLELAB’s detection and PyHFO. The overall discrepancy was defined as the sum of the number of new HFOs detected by RIPPLELAB (new-RIPPLELAB) and the number of new HFOs detected by PyHFO (new-PyHFO) divided by the total number of HFOs detected by RIPPLELAB. The discrepancies between PyHFO and RIPPLELAB of STE and MNI detectors were 10% and 14%, respectively.

[0162] As highlighted above, the implementation within PyHFO closely adheres to prevailing methods for data reading. Additionally, it provides a more accurate representation of the input parameters utilized in the construction of the bandpass filter. Consequently, PyHFO’s methodology exhibits greater implementation accuracy compared to most mainstream publicly released software.Table 2.3. Event comparison of differences between RIPPLELAB and PyHFO implementations in UCLA, Zurich, and Rodent dataset. The new-RIPPLELAB row represents the number of new events detected by RIPPLELAB, and the new-PyHFO row represents the number of new events detected by PyHFO for the specific experiment we conducted when the agreement of two events is defined as the 50% overlap of two events. Please note that the difference between 90% overlap and 50% overlap is minimal, amounting to no more than 0.2% (Noverlareading and bandpass filterwere from RIPPLELAB, but detection algorithm was from PyHFO; Exp2: data reading was from RIPPLELAB, but bandpass filter and detection algorithm were from PyHFO; Exp3: bandpass filter was from RIPPLELAB, but data reading and detection algorithm were from PyHFO; PyHFO: all data reading, bandpass filter and detection algorithm were from PyHFO.

[0163] HFO detector runtime comparison

[0164] Table 2.4 presented a runtime comparison between PyHFO and its Matlab-based counterpart (RIPPLELAB) across various hardware specifications on UCLA, Zurich, and Rodent datasets. To save computational resources, we reported only the runtime of RIPPLELAB on the Linux machine. Since different datasets had different numbers of recordings and lengths, we normalized the runtime to measure the detection speed: how many seconds the detector will take to process one minute of recording in one channel (120k data samples at a sampling frequency of 2,000 Hz). We also put the rough total time processing each dataset in parentheses (see Table 2.5 for detailed runtime in minutes). When comparing the detection speed of HFO detection using the same detection parameters, PyHFO significantly outperformed RIPPLELAB in both single-core (n=l) and multi-core (n>l) configurations, as detailed in 2.4.6 for hardware setup specifications. For the UCLA dataset, the detection speed of PyHFO could be at least 50 times faster (1.309 seconds / channel / minute for RIPPLELAB and 0.018 seconds / channel / minute for PyHFO on STE detector, and 19.82 seconds / channel / minute for RIPPLELAB and 0.292 seconds / channel / minute for PyHFO on MNI detector). PyHFO in the UCLA dataset (median of 94 channels per recording) can utilize multiprocessing (when n=32) better than in the Zurich dataset (median of 23 channels per recording) and Rodent dataset (median of 9.5 channels per recording). Nonetheless, PyHFO ran at least 15 times faster than RIPPLELAB in these three datasets. Even when operating with a single core, PyHFO still offers at least a six-times improvement in speed. Compared to STE, the MNI detector’s longer runtime is due to an iterative procedure within its computational pipeline. It took days for RIPPLELAB to detect MNI HFOs for a fairly large dataset (4 days for the UCLA dataset and 6.6 days for the Zurich dataset), which blocked feasibility for large- scale HFO analysis in the communit '. However, using PyHFO under a multi-core setting, the runtime could be significantly reduced to within six hours.Table 2.4. Comparative analysis of runtime (measured in runtime (seconds) / channel / recording minutes) in RIPPLELAB and PyHFO: Detection of all events from UCLA, Zurich, and Rodent dataset. We put the roughly total time of detection into parentheses. When n>l represents when PyHFO runs in a multi-core setup, n-jobs=32 for the Linux machine and n-jobs=8 for macOS and Windows machines. The best performance in each machine and dataset was highlighted in bold. Abbreviation: d: day(s), h: hour(s), m: minute(s), s: second(s).Table 2.5, Comparative analysis of runtime (measured in minutes) in RIPPLELAB and PyHFO: Detection of all events from UCLA, Zurich, and Rodent datasets. When n>l represents when PyHFO runs in a multi-core setup, n-jobs=32 for the Linux machine and n- jobs=8 for macOS and Windows machines. The best performance in each machine and dataset was highlighted in bold.

[0165] HFO event annotation

[0166] The HFO annotation was conducted on STE HFO events from the UCLA dataset (n = 12,494). Two experts (NHand SH) annotated HFO events into one of the three classes: artifact, HFO-with-spike (spkHFO), and HFO-without-spike (non-spkHFO). The inter-rater reliability of these two expert annotators was measured by the Cohen kappa score (kappa = 0.96 for labeling artifact, 0.85 for labeling HFO-with-spikes). The evaluation procedure was reported in Zhang. This annotation yielded 6,294 HFOs with spikes (spkHFO), 3,459 HFOs withoutspikes (non-spkHFO), and 2,741 artifacts. To ensure the DL models trained from the STE detector also generalize well in HFO events detected by the MNI detector, an expert (HN) annotated MNI HFOs into the artifact, spkHFO, and non-spkHFO from representative subjects (3 subjects, n = 758, including 416 artifacts, 312 spkHFOs, and 30 non-spkHFO) and the performance metric of the model against the expert annotation was reported.

[0167] Machine learning algorithm against expert labeling

[0168] In five-fold cross-validation, for STE HFOs (19 subjects, n = 12,494), the model achieved an accuracy of 98.6% and 89.1% for classifying artifacts and HFO with spikes, respectively, as shown in Table 2.6. This performance is almost the same as that in Zhang (artifacts: 98.8%, spkHFO 89. 1%) but with much lower MACs and runtime when we classified HFOs in CPUs in Table 2.7. More importantly, the trained model using STE HFOs could successfully classify the MNf HFOs, demonstrating the success of the data augmentation and generalization ability of the model, which enables these two deep learning models to be used in natural settings. The excellent performance across detectors also demonstrates the morphological similarity between the spkHFO in MNI and STE detectors.Table 2.6. Performance analysis using five-fold cross-validation: mean of accuracy, fl -core, recall, and precision on the test set in five-fold cross-validation with 95% confidence interval versus expert labeling.Table 2.7. Comparative analysis of computational costs across models: MACs and model run-time (in seconds) for GPU, Linux, mac OS, and Windows (CNN was the model design in Zhang; Dim. Optim. was the run-time after input dimension optimization; Pruning was the model runtime after dimension optimization and pruning). The best performance was highlighted in bold.

[0169] Neural network complexity comparison

[0170] In Table 2.7, we compare the Multiply Accumulate (MACs) on a single data sample as input and runtime of inference for 1000 data samples using GPU and CPUs between state- of-the-art and PyHFO. We reported the performance metrics of spk-HFO and artifact classifiers, respectively. By computing the MACs, the classifiers in PyHFO are more computationally efficient than the models proposed in Zhang, which provides theoretical support for later empirical experiments. Furthermore, even though both classifiers from Zhang and PyHFO run at comparable speeds in GPU, the artifact classifier in PyHFO runs at least 4 times faster than its counterpart, and PyHFO ’s spkHFO detector runs three times faster in CPUs. As another ablation study, we blindly pruned the published network using the same pruning and fine-tuning parameters but without input dimension optimization. Even though the performance of the pruned model was still comparable, MACs were still around 500M, which is much higher than our approach.

[0171] Discussion

[0172] Our work was implemented based on strong clinical motivation. Prior observational studies and a clinical trial have shown issues with time constraints in HFO analysis in clinical settings. The clinical use of HFO detection sites during epilepsy surgery planning requires a fast, reliable, and user-friendly application of HFO detection. It also needs to simulate human experts’ judgment to complement the entire process, including the detection and classification of HFOs. Our platform incorporated such capacities and has the ability to use multiple HFO detection methods and classifiers, including artifact rejection and HFOs with spikes vs. without spikes. Additionally, this system is portable, allowing any physician or researcher with a laptop to utilize DL-based algorithms in various settings, such as the epilepsy monitoring unit or the operating room. This capability has the potential to facilitate clinical trials. While developing our PyHFO application, we demonstrated that our HFO detection algorithm is comparable with another open-source work, RIPPLELAB. We comprehensively tested our HFO detection algorithms on three datasets: UCLA (grids / strips), Zurich (grids / strips and depth electrodes), and Rodent dataset. While we implemented the Python version of this HFO analysis application, EEG data reading and input format were deployed using the Python package. We followed the same EEG reading calibration as MNE, while in RIPPLELAB, the calibration was only done by voltage, without the offset adjustment. Furthermore, there were slight differences in the data filtering implementation. We chose to use the filter construction by Scipy as it can produce a more accurate frequency response. Our study reported minor differences in HFOdetection numbers between RIPPLELAB and PyHFO, and we concluded that those are based on differences in data reading and filtering implementation between Matlab and Python.

[0173] We combined multiple methods to decrease the runtime of the whole pipeline. We utilized the multiprocessing feature of Python, employed vectorization implementation in wavelet computation, optimized the neural network’s input size, and pruned the neural network architecture. We also developed a data augmentation strategy in the neural network training to improve the generalization ability of the model. We demonstrated that with the use of our application, the runtime was at least 15 times faster in STE detection and MNI detection compared to RIPPLELAB. We also achieved high performance in classifying artifacts and HFOs with spikes (98.6% and 89.1%, respectively, on five-fold cross-validation and 98.7% and 89.8%, respectively, on an independently annotated test set). Our study did not validate our detected HFOs against clinical outcomes, such as postoperative seizure outcomes. Rather, we aimed to establish the engineering validity of our detection algorithms using grid / strip, SEEG, and rat EEG data.Example 2

[0174] Methods

[0175] Patient Cohort: This was a multi-institutional retrospective cohort study. The inclusion criteria consisted of [a] simultaneous video-iEEG recording for epilepsy surgery between August 2016 and December 2023 at UCLA Mattel Children's Hospital, or between January 2007 and May 2018 at Detroit Medical Center, [b] iEEG sampling rate of at least 1,000Hz, [c] iEEG contained at least an artifact-free 20 min slow-wave sleep epoch at least two hours apart from clinical seizure events, and [d] known postoperative seizure outcomes over one year. The exclusion criteria included [a] undergoing hemispherectomy or hemispherectomy, and [b] the presence of massive brain malformations (such as megalencephaly and perisylvian polymicrogyria) or previous surgeries that make it difficult to identify brain anatomy during the iEEG study. The institutional review board at UCLA and Wayne State University have approved the protocol. We obtained written informed consent from patients or the guardians of pediatric patients.

[0176] Patient evaluation: All patients with medically refractory epilepsy referred during the study period underwent a standardized presurgical evaluation, which — at a minimum — consisted of inpatient video-EEG monitoring, high resolution (3.0 T) brain magnetic resonance imaging (MRf), and 18 fluoro-deoxy glucose positron emission tomography (FDG-PET), with MRf-PET co-registration. The margins and extent of resections were determined mainly basedon the SOZ, clinically defined as regions initially exhibiting sustained rhythmic waveforms at the onset of habitual seizures. In some cases, the seizure onset zones were incompletely resected to prevent an unacceptable neurological deficit.

[0177] iEEG recording: Macroelectrodes, including platinum grid electrodes (10 mm intercontact distance) and depth electrodes (platinum, 5 mm intercontact distance), were surgically implanted. The placement of intracranial electrodes was mainly guided by the results of scalp video-EEG recording and neuroimaging studies. All electrode plates were stitched to adjacent plates, the edge of the dura mater, or both, to minimize movement of subdural electrodes after placement. Regarding the SEEG placement, BrainLab Elements software was used for planning the electrodes to the intended targets using T1 -weighted sequences, and the trajectories were guided by a gadolinium-enhanced T1 -weighted MRI. iEEG recording was obtained using Nihon Kohden Systems (Irvine, California, USA) for both institutions. The sampling frequency was set at 1,000 Hz in Detroit and at 2,000 Hz in UCLA upon acquisition.

[0178] Acquisition of three-dimensional (3D) brain surface images: We obtained preoperative high-resolution 3D magnetization-prepared rapid acquisition with gradient echo (MPRAGE) T1 -weighted image of the entire head. Using the FreeSurfer scripts, we created the averaged surface image to which all electrode locations were spatially normalized. In cases where the software failed to detect the pial surface accurately due to insufficient cerebral myelination, we manually delineated the pial surface using the Control Point function. For patients younger than two, we used the Infant FreeSurfer software package to reconstruct the surface image. The averaged surface image functioned as the template for the analysis of anatomical location.

[0179] Anatomical labeling and determination of ROIs: For the dataset from UCLA, each implanted contact was labeled visually according to the Desikan-Killiany-Tourville atlas. The location of electrodes was directly defined within a Freesurfer-based 3D surface image using post-implant computed tomography (CT) images using Brainstorm software. For the dataset from Detroit, all implanted subdural contacts were coregistered with 3D surface images within the FreeSurfer with an FSaverage vertex label. We defined 34 ROIs for further analysis (Table 2.8). For the data harmonization between the two institutions, the FSaverage vertex of Detroit datasets was converted to MNI coordinates. Finally, these data were combined with UCLA patients, which were projected to the MNI normalized space under Brainstorm for the coregistration image (FIG. 9A).

[0180] iEEG data pre-processing: We used a customized common average reference for the grid and a bipolar montage for SEEG data. The EEG was first resampled to be the samesampling frequency of 1,000 Hz, and a band-reject filter was to reject the 60 Hz and its harmonics with a stopband of 2 Hz. iEEG channels not recording from grey matter (e.g., located outside of the brain) or otherwise deemed 'bad' (e.g., excessively noisy or artifactual) by the clinicians were discarded from the analysis. All EEG data pre-processing and analysis were performed using Python 3.9.1 (Python Software Foundation, Wilmington, DE).

[0181] Automated HFO detection: HFOs were detected by the PyHFO platform using both STE and MNI detectors to enhance the sensitivity of detection. For STE and MNI detectors, the detection parameters were set to default, except that the frequency band for filtering the signal was set to be 80-300 Hz to accommodate a minimum 1,000 Hz sampling frequency.

[0182] Overall training and inference method: The overall study flow is outlined in FIGS. 9A- 9D and FIGS. 10-12.

[0183] Morphology-based pathological HFO classification by deep generative model

[0184] Subject-wise k-fold Cross-validation'. To thoroughly test our method, we use a subjectwise five-fold cross-validation. In this setup, for each fold, we set aside 20% of the subjects as a test set (controlling sampling from UCLA grid / strip, UCLA SEEG, and Detroit grid / strip, ensuring that over five folds, every subject is tested once. All remaining data became the training set. Within the training set, we then randomly sampled ten subjects uniformly from three datasets for the validation set.

[0185] Feature Representation'. Each event detected by the automatic detector was represented by a time-frequency plot (Morlet wavelet transform), as the time-frequency plot is known for representing spatial -temporal information for each HFO event. The time-frequency plot spans ± 285 ms (centered in the middle of the HFO event) in time and 10-290 Hz in frequency. This plot was resized to a 64 x 64 resolution, and the value was normalized within the range of 0 to 1 before being sent to the VAE encoder.Table 2.8. Spatial distribution of intracranial electrode sampling.

[0186] VAE architecture: VAE emerges as a prominent deep learning mode. Distinct from traditional autoencoders, VAE introduces a probabilistic approach; the encoder encodes inputs as distributions. During the decoding process, the model stochastically samples from these distributions, enabling it not only to accurately reproduce known data but also requires the decoder to generate new similar data. The loss of training the VAE, reconstruction loss, and variational regularization, the reconstruction loss measures the distance from the input to the output, and the variational regularization ensures the distribution generated from the encoder follows the normal distribution. We've chosen ResNet as the backbone for both the encoder and decoder, drawing inspiration from its widespread application in the community of computer vision, as evidenced by sources, as well as the successful application of using ResNet in capturing morphological information within the time-frequency plot. The designated size forthe latent space is fixed at 8. As an ablation study, we set the latent dimension to 16 and demonstrated that it would produce comparable results but with more redundant latent space by visualization (FIGS. 13-14).

[0187] Self-supervised VAE training: Our VAE undergoes training over 80 epochs with the main goal of reconstructing the input image (FIG 9A). During the training, we used Adam optimizer with a learning rate 3e-4 and batch size of 512. Furthermore, the time-frequency plot is augmented by randomly flipping at the time axis to improve the model's generalization ability. To ensure an even representation of each subject during every epoch, we employ a stratified sampling method, capping the sample number at 2,500 per subject in each epoch of training. Moreover, we selected the perceptual loss as our reconstruction loss criterion to more effectively capture the morphological discrepancies between input and output images. We adopted the beta-VAE in the training to better disentangle the latent space; therefore the final loss function loss = perceptual loss + beta*KL divergence, where beta = 0.1. The model iteration (epoch) exhibiting the minimum validation loss was the one selected for subsequent utilization in the unsupervised classifier assembly and inference.

[0188] Unsupervised discovery of HFO clusters: We designed a hierarchical (two-stage) Gaussian Mixture Model (GMM) based clustering pipeline to discover morphologically different classed of HFOs (FIG. 9C). Once the VAE had been trained, we processed the entire training dataset through the VAE encoder, thereby extracting the associated latent codes (FIG. 9B). After this extraction, at the first stage, we employed the GMM clustering algorithms, setting the cluster number, k, to 2. We adopted stratified sampling to balance the contribution of information from each subject. We capped the number of events from each patient at 10,000. The cluster with high reconstruction loss was defined as the cluster of artifacts, as the artifact is diverse in morphology and the VAE cannot properly decode. Then, at the second stage, within the cluster that was classified as non-artificial events, we trained another GMM to cluster them further into two clusters. We capped the number of events from each patient at 2000. To allocate either a pathological (mpHFO) or physiological (non-mpHFO) label to each cluster, we employed a minimalistic use of clinical data. The cluster with a higher resection percentage in seizure-free patients after resection was deemed pathological. As part of an ablation study, we also explored different stratified sampling methods, which demonstrated similar performance in prediction surgical outcomes (FIG. 13).

[0189] HFO morphology inference pipeline: For predicting the HFO events from a new subject, time-frequency plots were first computed. Then, these time-frequency plot were sentinto the VAE encoder to extract the latent codes. The trained hierarchical GMMs predicted each latent codes and assigned class labels: mpHFO, non-mpHFO and mArtifact for each detected event.

[0190] mpHFO characterization

[0191] VAE-based HFO morphology inference pipeline:

[0192] Time-frequency plot characteristics of pathological and physiological HFOs: We investigated whether the time-frequency scalogram of mpHFOs identified by the VAE exhibited a distinction from that of non-mpHFOs (FIGS. 15D-15E). We used similar methodologies outlined in preceding studies.21Specifically, we executed a t-test across every pixel in the time-frequency scalogram comparing the predicted mpHFOs and non-mpHFOs. The foundational null hypothesis is that scalogram values associated with mpHFOs are higher than those linked with non-mpHFOs. We binarized the P-value less than 0.01 (if the P-value is less than 0.05, the value is 1, and vice versa). The characteristics map displays mpHFOs’ morphological distinction within frequency and time domains.

[0193] Latent space interpretability analysis:

[0194] Latent space 2D visualization: The encoder extracts latent vectors from the timefrequency plot of each event, and the latent code could contain essential morphological information that the decoder can reconstruct. In order to visually demonstrate the topology of the latent space trained from VAE better, we used t-distributed stochastic neighbor embedding (TSNE) to project the latent space acquired from the VAE into 2D space. We used the TSNE accelerated by GPU (cuml package), the parameters of fitting the TSNE using default parameters as the released package except for the n_components as 2 (FIG. 11. Model interpretation analysis workflow). We projected all of the latent code into the 2D space. For better visualization, we random sampled maximum 200 samples from each subject for visualization, and color code each project 2D points by different types of HFOs: mpHFO, non- mpHFO, mArtifact (FIG. 15A); spkHFO, non-spkHFO, artifact (FIG. 15B); different subcategories in subjects demographic information including data source (FIG. 16 A), gender (FIG. 16B), age (FIG. 16C), and different anatomical locations (FIG. 16D).

[0195] Statistical tests to evaluate the dependence of HFO morphology on the significance of different demographic subject-level information and anatomic locations of the HFO events variability:

[0196] We evaluated if HFO morphology is dependent on demographic and anatomical factors using statistical tests. Specifically, we assessed whether variables such as gender influencedHFO characteristics. We used logistic regression classifiers, training them with latent codes of HFOs labeled by gender. We selected three male and three female subjects per fold, sampling 100 HFO events from each to form a training set of 200 latent codes. We repeated this for five trials per fold, totaling 25 accuracy data points for each classifier. For a comparative analysis, we implemented a surrogate classifier with randomized gender labels. Using the same sampled latent codes, we shuffled and reassigned the gender labels randomly across the patients, thus maintaining the data integrity but altering the label associations. This classifier was also run for 25 trials per fold to ensure consistency in comparison. Additionally, we extended this methodology to investigate the impact of anatomical locations on the morphology of both putative physiological HFOs (from non-seizure regions in seizure-free patients) and putative pathological HFOs (from seizure onset zones, SOZ). Each anatomical location was uniquely coded, and we sampled 100 HFO events from each location for every patient in a fold. This resulted in 500 events per trial for training a multi-class logistic regression model. We conducted five trials per fold, generating 25 sets of accuracy measures and confusion matrices. We compared the true-label classifier's performance to that of the surrogate classifier using one-tailed t-tests. This statistical test determined if the accuracy of classifiers trained on true labels was significantly higher than those trained on randomized labels, thus indicating a dependency of HFO morphology on the tested variable.

[0197] Suppose HFO morphology depended on a particular variable (for example, the gender of the subject, FIG. 16C). In that case, the latent codes of the HFOs labeled with the corresponding subcategories (female or male) may be distinguishable. Consequently, a classifier designed to identify true labels from the latent codes of HFOs should achieve higher accuracy compared to a classifier trained with randomly assigned labels. The methodology can be illustrated using gender as an example. For each fold of the cross-validation, we randomly selected a sample of male (n=3) and female (n=3) patients. We randomly sampled 100 HFO events from male subjects and 100 HFO events from female subjects. This resulted in 200 latent codes for training a logistic regression classifier. For computing accuracy, we repeated the sampling process for the left-out subjects for that fold. This resulted in one accuracy datapoint. We repeated this sampling, training, and testing process five times (trial) for each fold. Thus, we have 25 trails in total (5 folds x 5 trails) for a classifier to determine gender based on actual gender labels. Next, we describe the design of a surrogate classifier based on randomizing the labels. For each fold and each trial, the same latent codes sampled for the true- label case are used. However, the gender label is randomized. That is, we randomly shuffle thesubjects and assign the label of male to the first three and female to the last three. Then, a surrogate logistic classifier is trained with the randomized gender labels. For testing, we again used the same data samples from the true-label case and randomized the gender assignment. This leads to another 25 accuracy samples for the surrogate classifier. A one-tailed t-test was employed to assess whether the true-label classifier's accuracy was significantly better than that of the surrogate.

[0198] We performed two such tests, one for the dependence of the morphology of HFOs from preserved regions in seizure-free patients (presumed physiological HFOs) on anatomical locations. The second one was for the dependence of the morphology of the HFOs from SOZ regions (presumed pathological HFOs) on anatomical locations (FIG. 17). Note that each channel in each subject was assigned a unique anatomical location. Given an HFO type, for each fold and each trial, the HFO events were sampled equally (n=100, in total 500) from each anatomical location across all patients. A multi-class logistic regression model is trained, and the confusion matrix was computed based on 500 samples acquired from the test subject under the same fold. Five such trials are done for each fold. Thus, for each HFO type, we obtained 25 confusion matrices, and the average of the confusion matrices and the accuracy were computed. Then, for the surrogate classifier, the randomization of labels of the 500 samples was done by shuffling samples on the channel level; for example, all samples from a particular channel that were initially labeled as "Frontal" would be reassigned to the same randomly selected label such as "Occipital". This led to 25 surrogate accuracy samples. A one-tailed t- test was employed to assess whether the true-label classifier's accuracy was significantly better than that of the surrogate.

[0199] Latent space disentanglement on latent dimension perturbation:

[0200] By manipulating specific dimensions of the latent space and observing the changes in the reconstructed images, we can infer the significance of each dimension in relation to the data characteristics it influences (FIG. 18, FIG. 11 : Latent dimension interpretation). This approach not only can aid in the interpretability of complex models but also provide insights into the underlying data structure that the VAE has learned to encode. The latent code was a onedimensional vector sampled from a distribution generated by the encoder using a timefrequency plot. The distribution of latent codes was plotted on each dimension and classification type (mpHFO and non-mpHFO) to identify the dimension that contributes the most to the decision boundary. To further understand the morphological characteristics contributed from each dimension, we manipulated each dimension of the latent code and thenobserved the resultant images generated from the decoder. The trend of the resultant images corresponding to the manipulation could reveal the morphological characteristic for that dimension. Specifically, for one latent code, we moved the value of a specific dimension to the 99th percentile and 1th percentile of all of the values within this dimension (we avoided extreme values to exclude outliers) while keeping the values in other dimensions unchanged. For visualization, we focused on the morphological impact on the two most representative latent codes, which were determined by averaging the latent vectors of all identified mpHFO and non-mpHFO events, respectively. To evaluate the effect of each dimension at a population level, we first randomly selected 500 high-confident mpHFO (model confidence > 0.99) and 500 high-confident non-mpHFO (model confidence < 0.01). We then manipulated the latent codes of 1000 randomly selected samples, irrespective of the event types (mpHFO / non- mpHFO). To quantitively evaluate the morphological change of resultant images, we designed various metrics based on simple image processing techniques on the time-frequency plot as a proxy to represent the morphological characteristic of the corresponding HFO events. To identify if an event was more pathological (determining pathological dimension), we calculated the ratio of power in the "hanging bell" template (FIG. 15E) to the overall power within the time-frequency plot. To calculate the peak frequency of an event (determining the peakfrequency dimension), we averaged the information across -90 ms to 90 ms (as the center of the event to be 0 ml) of the time-frequency plot to construct a frequency-power vector and found the frequency that had the max power within 80 Hz to 290 Hz. To measure the slow- wave power of an event (determining the slow-wave dimension), we computed the power in the 10-20 Hz range on the time-frequency plot to the total power. These analyses were visualized using box plots. Furthermore, we used the inference pipeline (encoder and GMM) to classify each generated event and visualized the model confidence through the violin plot.

[0201] Clinical correlation: Predicting surgical outcomes: Once the VAE and GMM classifiers are trained on the training set, this pipeline can make a prediction on the test set. In the following analysis, we use subjects whose number of detected events greater than 50 per recording (n=l 63). We evaluated the effectiveness of a discovered biomarker by predicting the surgical outcomes of the resected subjects. Specifically, features associated with each subject are extracted, and then a predictive model is trained using these subject-specific features to predict whether or not one subject will be seizure-free. In particular, we used resection ratio (No. resected events / No. overall detected events), which is the same as in many HFO studies and subjects' demographic features such as gender and age. Furthermore, we also consideredthe resection status of SOZ, which is a current clinical standard to guide epilepsy surgery. We applied two approaches to validate and predict the surgical outcomes. First, we evaluated the separation of post-surgical seizure-free subjects and post-surgical none-seizure-free subjects in the feature space and also evaluated the balance of the trained unsupervised models in the fivefold cross-validation, i.e, if a model trained in one fold is too conservative or aggressive on predicting the biomarkers. Therefore, we trained a single / multivariable logistic regression model to predict post-surgical outcomes on features generated from all subjects in all test sets among the five-fold cross-validation and compared the AUC of the logistic regression. Second, we evaluated the forward predictive power of our discovered biomarker. In each fold of fivefold cross-validation, we trained a random forest model from the training set and validation set, using subject-associated features as input of random forest and surgical status as labels. We used the trained random forest to predict all subjects underwent resection in the test set. Then, one random forest model was trained for each fold, and a set of performance metrics (accuracy, precision, recall, and Fl score) were reported; we reported the mean and SEM of the Fl score based on the unbalanced dataset.

[0202] Statistical analysis: The above-mentioned statistical calculations were carried out using Python (version 3.9.1; Python Software Foundation, USA). The deep neural network was developed using PyTorch (version 2.1.0; Facebook's Al Research lab). Quantitative measures are described by medians with interquartile or mean with standard deviations. Comparisons between groups were performed using chi-square for comparing two distributions and Student's t-test for quantitative measures (in means with standard deviations). All comparisons were two- sided, and significant results were considered at p < 0.05 unless stated otherwise. Specific statistical tests performed for each experiment were described in each section. Machine learning model performance was evaluated using accuracy ([TP + TN] / [TP+TN+FP+FN]), recall (TP / [TP+FN]), precision (TP / [TP+FP]), and F-1 score (2 / [l / recall + 1 / precision]).

[0203] Results

[0204] Building an explainable unsupervised DL algorithm for HFO analysis:

[0205] We studied 185 patients (91 females) from two centers who met the eligibility criteria (Table 2.1). The median age at surgery was 13 years (range: 2-44 years). A total of 18,265 artifact-free electrode sites (median: 106 per patient; range: 29-152) within the regions of interest (ROIs) were available for analysis (FIG. 9A). There were 1,670 electrode sites sampled within the SOZ (mean: 9.42 per patient; frontal: 437; temporal: 464; parietal: 393; occipital: 135; limbic: 241), and 7,732 sites sampled within non-epileptogenic brain regions (frontal:2943; temporal: 1806; parietal: 1853; occipital: 635; limbic: 495) (Table 2.8). The median duration of analyzed EEG data for the UCLA grid / strip dataset was 91.5 min [IQR: 90.6-94.8 min], and for SEEG, it was 90.3 minutes [IQR: 87.1-96.8 min]. While for the Detroit dataset, the median analyzed EEG recording duration was 5.3 minutes [IQR: 5. 1-5.7 min]. In total, 686,410 putative HFOs were detected from all the datasets by the PyHFO platform using both STE and MNI detectors.52The median rate of HFOs (number of detections / min) within SOZ contacts was 3.11 (range: 0.73-7. 17), and the median rate of HFOs within the non-epileptogenic contacts was 1.08 (range: 0.63-2. 15) across the ROIs. All the detected HFOs were used for the unsupervised VAE model training (FIG. 9B). Using the trained VAE model, such HFOs were first classified into artifacts (mArtifacts) and non-artifactual HFOs using GMM based on reconstruction loss (FIG. 9C). All the non-artifactual HFOs were further classified into two clusters. The cluster with a higher resection percentage in seizure-free patients after resection was deemed pathological. HFOs coming from the pathological cluster were defined as mpHFOs. All the detected HFOs were also classified into artifacts, HFOs with spikes (spkHFOs) and HFOs without spikes (non-spkHFOs) using our previous classification algorithm for the subsequent correlational analysis, including interpretability analysis (FIG. 9D). The detailed workflow of the study is described in FIGS. 10-12. There were 163 patients who underwent resective surgery, and 110 patients (67.5 %) achieved seizure freedom. Of the patients who had resection, pathology results were as follows: focal cortical dysplasia (FCD) (41.1%), hippocampal sclerosis (HS) (6.7%), tumor (19.0%), and others (33.1%).

[0206] Characterization of pathological HFOs based on the unsupervised VAE algorithm:

[0207] With five-fold subjective- wise cross-validation, the VAE classification was aligned with conventional knowledge that spkHFOs are more likely pathological: 93.56% of mArtifact were predicted as artifactual HFOs, 80.1% of the mpHFOs were predicted as spkHFOs and 84.87 % of the non-mpHFOs were predicted as non-spkHFOs (p < 0.001). Furthermore, we visualized these findings by projecting the latent codes of each event onto a two-dimensional plane using t-distributed Stochastic Neighbor Embedding (t-SNE). (FIGS. 15A-15B). Such clustering results were consistent throughout the five folds (FIG. 19). Moreover, the mpHFO rate (n_event / min / channel) was significantly higher than the discovered non-mpHFO rate (n_event / min / channel) in SOZ channels across three datasets (one-tailed t-test, p-value < 0.01 for UCLA grid / strip dataset, p-value < 0.01 for Detroit grid / strip dataset, and p-value < 0.01 for UCLA SEEG dataset) and did not exhibit significant difference within the non-soz channels (two-sided t-test, p-value = 0.1 for Detroit grid / strip dataset, p-value = 0.49 for UCLA SEEGdataset, and p-value = 0.38 for UCLA grid / strip dataset) (FIG. 15C). We then investigated the neurophysiological characteristics of mpHFO and non-mpHFO directly from the data. The morphological analysis of the time-frequency plot demonstrated that mpHFOs had higher values throughout the HFO band (> 80 Hz), around the center point (0 ms, where HFOs were detected) than non-mpHFOs (FIG. 15D-15E). Furthermore, there were statistically higher values of mpHFOs at the low-frequency (non-HFO) band throughout the time window compared to non-mpHFOs. These bands together lead to a "hanging bell" template in the timefrequency plot (FIG. 15E). Such a template showed consistency regardless of the variables, including gender, the origin of the dataset, pathology, and age categories (FIG. 20). mpHFOs also showed distinct characteristics in the frequency domain with FFT analysis. The mpHFOs demonstrated relatively high amplitude within the low frequency band, while non-mpHFOs, likely representing physiological HFOs, exhibited discrete peaks with the HFO range, maximal at 22.8 Hz (FIG. 15F-15G). To understand the low-frequency band component of the "hanging bell" template in mpHFOs, we filtered mpHFOs with a 10-80 Hz bandpass to evaluate the corresponding time-series data. We found a spike-like EEG signal after this procedure (FIG. 15H-15J) but none in non-mpHFOs (FIG. I5K.-15M).

[0208] HFO morphology analysis based on dataset origin, sex, age, and pathology:

[0209] We investigated whether HFO morphology was affected by variables including recording sites / type (UCLA grid / strip, UCLA SEEG, Detroit grid / strip), gender (female vs. male), age (0-5, 6-10, 11-15, 16-20, and 21+), and pathology (HS, FCD, Tumor, others). We first visualized the distribution of the latent space by projecting latent codes from each subject into 2D space using t-SNE and color-coded with different subcategories within that specific demographic variable (FIG 16A, 16C, 16E, 16G). The projected codes were uniformly randomly sampled 300 from each subject to achieve easier visualization. Then, a statistical test was conducted by comparing the accuracy of a logistic regression model trained from actual data with the accuracy of a surrogate logistic regression model trained with label-permuted data. HFO morphology' based on the dataset sources could not be differentiated from the classifier (p-value = 0.13, FIG. 16B). Similarly, HFO morphology, based on demographic information such as gender, age group, and pathology, could not be differentiated (Gender: p- value = 0.44; Age Group: p-value = 0.10; Pathology: p-value = 0.83) (FIG. 16D, 16F, 16H).

[0210] HFO morphology analysis based on ROIs (anatomical location):

[0211] We then investigated whether HFOs generated from different anatomic locations exhibit distinguishable morphologies (frontal, temporal, parietal, occipital, and limbic).Initially, we generated a plot of the latent space for HFOs derived from the preserved regions in patients who remained seizure-free after resection. The assumption was that such HFOs were likely physiological. We found a cluster of HFOs corresponding to the occipital region, which were distinguishable from other areas (FIG. 17A-17B). The confusion matrix showed that such occipital HFOs could be successfully distinguished by the classifier, which achieved an average of 62% accuracy (FIG. 17C). To further characterize such distinguishable features, we plotted the summation of the time-frequency plot of each HFO event in different anatomical locations (FIG. 17D). The HFOs derived from the occipital area showed distinct peak frequency and power ratio compared to HFOs generated from other anatomical areas (FIG. 21). Nevertheless, HFOs generated from SOZ did not demonstrate distinctive morphological differences across the anatomical locations using the same method (FIGS. 17E-17H).

[0212] Disentanglement of the latent space to establish neurophysiological characteristics of pathological HFOs:

[0213] We first plotted each dimension of the latent space and colored them with mpHFO and non-mpHFO for a specific fold (FIG. 11 : Latent dimension interpretation); the mpHFO and non-mpHFO exhibit different distributions in some of the dimensions but do not in the others (FIG. 22). We then used the decoder of the VAE to understand the meaning of each latent dimension by perturbing each dimension in the latent space and visualizing the reconstructed image. Specifically, in each fold, we first took the mean of the cluster of the mpHFO and non- mpHFO and then perturbed a single dimension from the smallest value (1 percentile) of that dimension to the largest value of that dimension (99 percentile). We found a dimension that demonstrates the separation of mpHFO and non-mpHFO. The morphological transition of the generated images showed more power within the "hanging bell" template region for both mpHFO and non-mpHFO when we traversed the value of that dimension (FIG. 18, panels a- b). At a population level, we randomly sampled actual latent code from both mpHFO and non- mpHFO and traversed that specific dimension of each latent code from the smallest (1% tile) to the largest (99% tile); the decoded images showed an increase of the power within the template region (FIG. 18, panel c). We fed the decoded time-frequency plots through our inference pipeline and plotted the distribution of the model confidence (FIG. 18, panel d), demonstrating the increasing model confidence toward mpHFO as we increased the value of that dimension in the latent space. Furthermore, we also captured dimensions that, while not separating mpHFOs and non-mpHFOs, still manifested interpretability neurophysiological properties. Within the same fold, a dimension was pinpointed correlating to the HFO's peakfrequency (FIG. 18, panels e-f), and subsequent population-level analysis indicated that an ascent in this dimension's value inversely correlated with the peak frequency in the timefrequency plots (FIG. 18, panel g) while the model confidence did not change towards either the mpHFO and non-mpHFO (FIG. 18, panel h). Similarly, the slow-wave dimension was revealed (FIG. 18, panels i-j), where an increase in dimension value directly paralleled an increase in power within the slow- wave band (10-20 Hz) in the decoded images (FIG. 18, panel k), but an increase in power within that dimension did not affect the model confidence towards either the mpHFO and non-mpHFO (FIG. 18, panel 1).

[0214] Prediction of postoperative seizure outcomes using the classified HFOs:

[0215] Based on the hypothesis that mpHFOs are pathological and their resection status can forecast postoperative seizure freedom, we investigated two prediction models. We first applied a logistic regression model with five-fold cross-validation. Notably, utilizing the resection ratio of the mpHFOs (AUC = 0.65) demonstrated superior separation capability compared to using the resection ratio of unclassified HFOs (AUC = 0.54) and traditional expert-driven classification, spkHFO (AUC = 0.60) (FIG. 14A). Further, a multivariable logistic regression model, incorporating subjects' demographic data as the baseline (age as a continuous variable and gender as a categorical variable) and the status of SOZ resection during surgery, was fitted. The mpHFO resection ratio again showed enhanced classification performance over the spkHFO (AUC = 0.71 vs. 0.68) (FIG. 14B). Additionally, we explored the feasibility of constructing a non-linear predictive model using a random forest classifier (FIG. 14C). The random forest aimed to forecast post-surgical seizure freedom in the test set based on subject-specific features identified in the training set through a subject- wise five-fold cross-validation, ensuring that our classifier was evaluated under consistent and rigorous conditions (FIG. 12). Our findings revealed that the random forest model, trained exclusively with the resection ratio of mpHFOs, exhibited better predictive performance (Fl = 0.73) compared to using unclassified HFOs (Fl = 0.66, p-value < 0.01) and spkHFO (Fl = 0.65, p- value < 0.01). The model's predictive power improved when incorporating more features, such as subjects' demographic information and SOZ resection status. Specifically, the combination of demographic data and mpHFO resection ratio (Fl = 0.76) outperformed traditional clinical predictions based on demographic information and SOZ resection status (Fl = 0.72, p-value < 0.01). Furthermore, a comprehensive model including all features, including demographic data, SOZ resection, and mpHFO resection ratio (Fl = 0.80) demonstrated superior predictive powerover the traditional expert-driven HFO classification including demographic data, SOZ resection, and spkHFO resection ratio (Fl = 0.78, p-value < 0.01).

[0216] Discussion

[0217] Leveraging the innovative self-supervised deep learning approach, we aimed to characterize pathological HFOs using a cohort of 185 epilepsy patients who underwent intracranial monitoring with grid or stereotactic EEG electrodes across 18,323 brain contacts from 34 ROIs. We analyzed 686,410 HFOs obtained via fully automated detectors. The novel VAE method and extensive data allowed us to control for variables including dataset origin, sex, age, epilepsy pathology, and anatomical location. By analyzing the latent space of VAEs, we identified morphological features of pathological HFOs, notably high signal intensity within the HFO band at detection, extending across beta and gamma frequency bands. These features remained consistent across all studied variables. Incorporating the resection status of the identified pathological HFOs into prediction models markedly enhanced the accuracy of postoperative seizure freedom forecasts, outperforming traditional approaches that utilize the resection status of SOZ.

[0218] Our results are compelling, as the morphologically defined pathological HFOs identified by our deep generative model exhibited features commonly associated with expert- acknowledged pathological HFOs, such as spikes. While "spike" annotations can be performed by human experts or computationally, there is no standardized method for determining what constitutes pathological spikes. Our deep generative model objectively identified distinctive features of pathological HFOs with a peak frequency of 23 Hz. Such non-HFO band EEG signals associated with pathological HFOs seemed to emerge as spikes in time series data. Rather than arbitrarily define spikes to analyze HFOs with spikes, this data-driven approach to characterize the "spike" component associated with pathological HFOs seemed more logical and unbiased. Another evidence of the pathological nature of mpHFOs was that they were primarily localized around SOZ, although the model was trained without a SOZ label.

[0219] Our findings provide clinically significant insights; the resection ratio of mpHFOs outperformed traditional unclassified HFOs and spkHFOs in predicting postoperative seizure freedom. This supports the idea that mpHFOs are more indicative of pathological activity than unclassified HFOs and spkHFOs. The potential clinical utility is evident as the resection ratio of mpHFOs exceeded the current clinical standards regarding the removal status of SOZ in predicting postoperative outcomes. mpHFOs can be analyzed from a short interictal EEG data obtained on the first night of recording, offering the potential to reduce the duration of EEGmonitoring, hospital stays, and associated costs for the patient. Furthermore, combining the resection status of SOZ with the proportion of mpHFO resection further enhanced the prediction, suggesting an additive effect when mpHFOs are combined with the current clinical standard.

[0220] Another notable discovery was that non-pathological HFOs originating from the occipital lobe displayed distinct morphological features compared to HFOs from other brain regions. Physiological HFOs were reportedly abundant in the occipital lobe and also showed distinct coupling with slow waves. This study added essential findings in the literature to establish the unique morphology of HFOs originating in the occipital lobe.

[0221] Until now, the field of HFO research has encountered significant challenges. HFOs can be generated in healthy brain regions, adding complexity to their interpretation. Despite promising results from numerous retrospective studies, some studies, including a prospective study and a clinical trial failed to demonstrate the utility of HFOs in improving postoperative seizure outcomes. The major challenge in HFO research is distinguishing between pathological HFOs that need to be targeted for resection, and physiological HFOs which can be spared during resection. Attempts to differentiate these include identifying HFOs associated with spike-wave discharges, which are generally considered indicative of pathology. However, there is no uniform method to define HFOs with spikes, and experts often exhibit inconsistencies and poor inter-rater reliability when using visual annotations. Moreover, simple analyses of HFO features such as frequency, amplitude, and duration have been ineffective in differentiating pathological from physiological HFOs. While fast ripples (250-500 Hz) may more accurately delineate epileptogenic zones than ripples (80-250 Hz), they are detected less frequently. Strategies like adjusting the detection rate of HFOs to match region-specific normative values can help estimate the pathology of corresponding brain regions, it cannot reliably classify each HFO event as pathological or physiological. This uncertainty complicates clinical decision-making, especially when treatment decisions hinge on a limited number of detected HFOs.

[0222] In this study, we adopted a fundamentally different approach. Instead of pre-selecting potential pathological HFO features, we implemented a data-driven, self-supervised strategy using VAE to define pathological HFOs. In the field of HFOs, prior studies have explored the use of various ty pes of neural networks, including convolutional neuronal network, long shortterm memory, and transformer, to classify HFOs into pathological and physiological based on human-annotated labels. This reliance limits scalability' on large datasets and is constrained bythe subjectivity inherent in expert labeling. Our method is contingent on the premise that pathological HFOs possess distinct morphological characteristics that are separable from their physiological counterparts. This hypothesis is biologically plausible, given that HFOs with similar frequency ranges emerge via markedly different generative mechanisms: Pathological HFOs arise from abnormal synchronous burst firing, whereas physiological HFOs result from inhibitory synchronous postsynaptic potentials. These morphological differences, potentially too subtle for detection with traditional signal processing methods using clinical macroelectrodes, may be effectively discerned through our innovative approach. VAEs are a DL framework capable of efficiently learning low-dimensional representations of unlabeled data. Without providing any labels associated with the input, an unsupervised classification can be achieved by reconstructing the original high-dimensional data from their latent codes. It excels at uncovering clusters in the latent space, with each cluster representing a distinct mechanism underlying the data generation. Consequently, clusters automatically formed in the latent space of a VAE have the potential to capture various morphological features and classify artifacts and HFOs with or without spikes, as demonstrated in our study. The diverse HFO morphologies generated by these clusters possessed distinct traits, allowing for interpretation through direct visualization and interpolations in the latent space, which enabled us to interpret what the algorithm has learned. The algorithm provided intuitive and interpretable latent spaces, such as a pathological template domain (the "hanging bell" template on a timefrequency plot), peak frequency domain, and slow-wave domain. VAEs have been applied in various EEG research studies, but this marks the first application of VAEs to HFO analysis in a large-scale dataset with robust neurophysiological and clinical correlation.Example 3

[0223] 3.1 Materials

[0224] 3.1.1 Patient Cohort

[0225] This study employed a retrospective cohort design. We consecutively recruited 159 children (under the age of 30) who were suffering from medically refractory epilepsy — characterized by monthly or more frequent seizures and the inefficacy of more than three antiseizure medications — from the Pediatric Epilepsy Program at UCLA and Wayne State. These children had intracranial electrodes implanted for epilepsy surgery planning and subsequent cortical resection between August 2016 and August 2018.

[0226] 3.1.2 Standard protocol approvals, registrations, and patient consent

[0227] The UCLA’s institutional review board approved the involvement of human subjects and dispensed with the necessity for written informed consent. All the tests were regarded as clinically relevant for patient care, and all the EEG data used in this retrospective study was anonymized before it was extracted and analyzed. This study was not classified as a clinical trial and was not registered in any public registry.

[0228] 3.1.3 Patient evaluation

[0229] All children referred during the study period, who were suffering from medically refractory epilepsy, underwent a systematic pre-surgery evaluation. This included inpatient video-EEG monitoring, high-resolution brain magnetic resonance imaging (MRI) of 3.0 Tesla, and 18 fluorodeoxyglucose positron emission tomography (FDG-PET), coupled with MRI- PET co- registration. The boundaries and scope of the resections were primarily based on the seizure onset zone (SOZ), which is clinically identified as the areas initially showing continuous rhythmic waveforms at the beginning of typical seizures. In certain instances, to avoid causing unacceptable neurological damage, the seizure onset zones were not fully resected.

[0230] 3.1.4 Subdural electrode placement

[0231] Surgical implantation was performed for macroelectrodes, which included platinum grid electrodes with an intercontact distance of 10 mm, and depth electrodes, made of platinum with a 5 mm intercontact distance. The total number of electrode contacts for each patient ranged from 40 to 128, with a median of 96 contacts. The positioning of these intracranial electrodes was primarily guided by the results of scalp video-EEG recording and neuroimaging studies. All electrode plates were sewn to neighboring plates, the dura mater’s edge, or both, in order to minimize the movement of subdural electrodes post-placement.

[0232] 3.1.5 Acquisition of three-dimensional (3D) brain surface images

[0233] We collected high-resolution 3D magnetization-prepared rapid acquisition with gradient echo (MPRAGE) T1 -weighted images of the whole head preoperatively. Utilizing post-implant computed tomography (CT) images, we created a FreeSurfer-based 3D surface image with the locations of the electrodes directly identified on the brain surface [Zhang et al., 2022b], Fur- thermore, digital photographs were captured intraoperatively before the dura was closed to improve the spatial precision of electrode placement on the 3D brain surface. When the area was re-exposed for resective surgery , we visually verified that the electrodes had remained in place, matching the digital photo taken during the electrode implantation surgery.

[0234] 3.1.6 Intracranial EEG (iEEG) recordings

[0235] The intracranial EEG (iEEG) data was recorded with a digital sampling frequency of 2,000 Hz using the Nihon Kohden Systems (Neurofax 1100A, Irvine, California, USA). For each participant, a consistent 90-minute EEG segment was chosen, ensuring it was at least two hours before or after any seizures, prior to tapering of anti-seizure medications, and before cortical stimulation mapping. This mapping usually took place two days post-implant. All the iEEG data used in the study were extracted from the clinical EEG recording.

[0236] 3.2 Dynamic directed weighted intenctal networks

[0237] In order to capture the temporal dynamics of the interictal interactions amongst the channels, we split the interictal iEEG data into one second long non-overlapping data segments. For each iEEG segment, we construct a frequency domain synchronization network. Assuming the interictal segment is T seconds long, we will have T such synchronization networks (FIG. 23). We describe the construction of one such network in the following section.

[0238] 3.2.1 Frequency domain synchronization network

[0239] We model the interictal interactions amongst the iEEG channels in the frequency domain. FIG. 24 shows a flow diagram showing the steps involved in constructing the frequency domain synchronization networks: i) select a specific frequency band to filter the intracranial EEG(iEEG) segment and transform the filtered signal into if s frequency domain representation and ii) subsequently, utilize the frequency domain representation to build a synchronization network to capture the interictal interactions amongst the channels. In the following sections, we describe each step in detail.

[0240] 3.2.1. 1 Frequency bands of interest in Epilepsy

[0241] The analysis of frequency bands within intracranial EEG (iEEG) data has provided valuable insights into the pathophysiology of epilepsy, aiding in diagnosis, treatment, and research endeavors. For example, increased delta power may indicate abnormal neuronal synchronization, especially during epileptic seizures onginating from deep brain structures like the thalamus while alterations in alpha rhythms, such as alpha power suppression, have been observed in the preictal and ictal states, suggesting their involvement in seizure generation and propagation. Similarly, ripple band and fast ripple band activity has been hypothesized to be a biomarker of EZ. Although, the framework that we are presenting is general and can be used for any frequency bands of interest we will be studying synchronization in High Gamma, Ripple and Fast Ripple band (FIG. 25).

[0242] 3.2.1.2 Network construction

[0243] The network consist of N nodes where N is the number of channels in iEEG. For each channel, we compute it’s power and phase spectrum using Fast Fourier Transform (FFT). Since the goal of the network is to capture synchronous activity between channels, so we capture the synchronous activity between channels i and j through phase coupling coefficient, whereand K is the number of discrete frequency bins in the frequency band of interest. However, we want to capture interactions during episodes of excessive synchronization so we scale the phase coupling coefficient by the total combined power in channels i and J, to get the synchronization coefficient between channels i and J

[0244] Hence, high total power along with large phase coupling coefficient will lead to large synchronization coefficient between nodes i and J. Having computed the synchronization coefficient, the next step is to map it into edge weights. We want the mapping to satisfy the following two criterion:

[0245] Edge weights should be bounded between 0 and 1 as it will facilitate network property computations. Edge weights should be a monotonically increasing function of the absolute value of the synchronization coefficient because high synchronization coefficient should map to a large edge weight.

[0246] The following mapping from synchronization coefficient to edge weight satisfies both the criterion

[0247] To emphasize substantial synchronization, we retain only the top 10 percent (ranked in terms of edge weights) of edges and discard the rest by setting their edge weights to zero. The final step is to assign directionality to the edges. We use the total phase method to assign the directionality (FIG. 26). The total phase in node j is given bywhere K is the number of discrete frequency bins in the given frequency band. Analogously, the total phase in node i is given bywhere K is the number of discrete frequency bins in the given frequency band. If there is a greater total phase in node J than node z, that isthen the weighted edge points from node i to node j (z — > j). If there is a greater total phase in node i than node 7, that isthen the weighted edge points from node j to node i (J z).

[0248] Seizure onset zone identification using dynamic interictal network properties

[0249] Recall, that our goal is to model the interictal interactions amongst the iEEG / SEEG channels using a dynamic synchronization network and leverage the temporal dynamics of the network properties, both local and global, to accurately identify the EZ (SOZ U PSOZ). Specifically, we want to build a classifier that takes in a node feature vector, constructed from the interictal network property dynamics, and outputs the probability of the node being a part of the seizure onset zone (SOZ) (FIG. 27).

[0250] In the next section, we explain how we construct the node feature vector.

[0251] 3.3.1 Node feature vector

[0252] We can observe from the interictal network snapshots shown in FIGS. 28-30 that they enter into a state of large scale synchronization where a star like network evolves with either a SOZ node or a resected Non-SOZ node being the hub. Since epileptic seizures are believed to be caused by large scale synchronization effects, so these hypersynchronous states and the nodes driving these states are of particular importance to us for SOZ identification.

[0253] A network feature which has been used to study synchronization in networks [Bertrand and Moonen, 2013] is Fiedler eigenvalue which is the second smallest eigenvalue of the Laplacian matrix L, whereL = D - A,D is the degree matrix and A is the adjacency matrix. In our context of interictal networks, we also observe a sharp rise in Fiedler value during episodes of large scale synchronization as shown in FIG. 31.

[0254] The most basic node feature to consider for capturing hub like phenomenon is in degree and out degree. Also, we want to put particular emphasis on the node degrees during these episodes of large scale synchronization so we also add the scaled degrees as node features, in degree and out degree scaled by fiedler eigenvalue. This scaling will amplify the node degrees during episodes of such large scale synchronization.

[0255] Another feature to consider at a node level is Eigencentrality because it has been shown in epileptic network research that eigencentrality of a node is correlated to the rate of pathological events generated by that node during seizure episodes. It measures a node’s importance not just by the number of connections it has, but also by the importance of the nodes it is connected to. Nodes with high centrality are not only well connected but are also connected to other important nodes in the network. It is computed by solving the set of consistency equations,where Ci and Cj are the eigencentralities of nodes i and j respectively. In matrix form, the equation iswhere c is the eigenvector centrality vector and A is the adjacency matrix. By Perron-Frobenius theorem, c is well-defined for any connected network and / . is the largest eigenvalue of ATwith corresponding eigenvector c.

[0256] In summary, for one interictal network, each node in the network is represented by a 5 dimensional feature vector

[0257] Since, we have a sequence of interictal networks so each node in the network is represented by a 5 dimensional time series representing the time evolution of this node across the feature space. Therefore, we need to summarize the 5 dimensional time series to get a feature vector for a node. To gain some insight into what type of summary statistics we should look into, we have plotted the time evolution of out degree (FIG. 32) and eigencentrality of a node (FIG. 33) in the network. From the plots, it can be observed that the node features occasionally enter into an outlier zone (marked in red). Therefore, we need higher order statistics, in addition to mean and variance, to capture this kind of bursty phenomenon of node features in time. The two higher order summary statistics that we use are kurtosis and skewness. Kurtosis quantifies the peakedness or flatness of a probability distribution compared to a normal distribution. High kurtosis indicates that the distribution has heavy tails (presence of outliers) and a sharp peak, while low kurtosis indicates a flatter distribution with lighter tails (absence of outliers). Skewness measures the asymmetry of the distribution : positive skewness indicates that the right tail of the distribution is longer or fatter than the left tail (presence of samples with extremely high values), while negative skewness indicates the opposite (presence of samples with extremely low values).

[0258] Putting everything together, each node is represented by a 5 dimensional feature vector that is evolving over time. We summarize the time evolution of each feature through 4 sufficient statistics thus giving a 20 dimensional feature vector per node (FIG. 34).

[0259] As mentioned earlier, we study synchronization effects in three frequency bands: high gamma, ripple and fast ripple. Each frequency band results in a sequence of synchronization networks which in turn gives a 20 dimensional node feature. Concatenating the features across the 3 frequency bands we get a 60 dimensional feature vector per node (FIG. 35).

[0260] 3.3.2 Node labeling

[0261] We have described the node feature vector construction in the previous section. The only step left is to assign node labels before we can train and test our SOZ classification model. We assign a label of 1 to all the SOZ nodes. For the rest of the nodes we can assign a label of 0 but that will make the training set noisy for the following reasons. In seizure free patients there might be channels in the resected zone that might have the capacity to generate seizures. From FIG. 30 we can see one such node (Resected Non-SOZ) being the hub of a hypersynchronous state and most likely it is PSOZ. Therefore, we drop resected Non-SOZ nodes from training. In non-seizure free patients, there might be channels in both the resected and preserved zone that might have the capacity to generate seizures since the surgical outcome was poor. Therefore, we drop non seizure free patients from training. Putting everything together, we only consider seizure free patients for training and for those patients we label SOZ nodes as 1 and Preserved Non-SOZ nodes as 0 (FIG. 36).

[0262] 3.3.3 Training and Testing corpus for SOZ identification

[0263] We constructed the training corpus by selecting seizure free patients with long interictal recording. Hence, 17 UCLA seizure free patients, each with an interictal recording of duration 90 minutes formed the training corpus and the remaining 142 patients formed the testing corpus. We trained one random forest SOZ classifier model on the channels of 17 patients and used it to infer the status of the channels of the remaining 142 patients (Table 3.1). However, we also want to infer the status of the channels of 17 patients. A patient-wise elimination training strategy was adopted to validate the findings of the SOZ classifier for these 17 patients. In this training strategy, we trained 17 models in total: one random forest model each for the 17 seizure free patients (FIG. 37). Therefore, we are using an independent cohort of channels to validate the findings of the detectors, as there is no overlap between the training and test set.

[0264] Random Forest is a machine learning algorithm that constructs a multitude of decision trees at training time and outputs the mode of the classes of the individual trees. Each tree is trained on a random subset of the training data and also uses a random subset of features for splitting at each node, thereby reducing overfitting and increasing robustness and accuracy.Table 3.1 Training and Testing Coprus: 17 UCLA seizure free patients in the training corpus and the remaining 142 patients in the testing corpus resulting in one random forest model

[0265] 3.3.4 Results

[0266] Recall, from our SOZ classifier setup shown in FIG. 27, that it outputs a probability of the node being part of SOZ. We can threshold the probability scores to get the predicted node labels and compare them with the ground truth labels. In Table 3.2, we tabulate the performance metrics when the trained SOZ classifier was used to predict the status of all the channels except the resected Non-SOZ ones. Since the ground truth label distribution is highly biased with only around 10% of the channels having a ground truth label of SOZ so recall is a more important metric than accuracy. From Table 3.2, we can observe that the classifier is able to correctly identify most of the SOZ channels in the test set shown by the high recall of 98%. The few false negatives (i.e. ground truth SOZ channels predicted as Non-SOZ), comes from patients where the SOZ labeling by doctors was generous. For example, most of the false negatives came from patients DL and GM where each of them had more than 30% of the channels labeled as SOZ.

[0267] When the same trained classifier was used to predict the status of the resected Non- SOZ channels, then around 18% of them got predicted as SOZ with a high probability. This result is along expected lines because from FIG. 30 we know that these types of channels also emerge as hubs of hypersynchronous states like the SOZ channels.Table 3.2 SOZ classification results: The SOZ classiifer outputs a probability score and we use a threshold of 0.5 to make predictions: P (z 6 SOZ) > 0.5 i E SOZ. The classifier is able to correctly identify most of the SOZ channels in the test set shown by the high recall of 98%. Few false negatives come from patients where the SOZ labeling by doctors were very generous.

[0268] 3.3.5 Discussion

[0269] We have modeled the interictal interactions amongst the iEEG channels using a sequence of synchronization networks and leveraged the temporal dynamics of the network properties to accurately identify the channels that drove all the seizure episodes (SOZ) duringthe 2 week long hospital monitoring with over 98% accuracy. One of the hallmarks of the constructed networks is that they occasionally transition into a state of hypersynchrony with SOZ nodes being the hub of these hypersynchronous states and our model is able to identify them and exploit them for accurate SOZ identification. Therefore, we have built a model that can identify the drivers of almost all the seizure episodes that happened during the 2 week long hospital monitoring from interictal segments with duration as short as 5 minutes. This result can significantly shorten the invasive hospital monitoring time, which is fraught with risk of infection.

[0270] The same model also identifies another set of channels (around 18%) of having the capacity to generate seizures and we hypothesize that these channels might be PSOZ.

[0271] 3.4 Predicting post-surgical outcomes using interictal SOZ classification and network dynamics

[0272] In the last section, we discovered a set of channels that got predicted as SOZ with a high probability but did not drive seizure episodes during the 2 week long hospital monitoring. These channels might constitute the unobserved / inactive part of EZ known as Potential Seizure Onset Zone (PSOZ). Since there is no diagnostic modality to know whether the entire EZ (SOZ U PSOZ) has been identified and resected, so the only way to know is through surgical outcome evaluation: if the patient becomes seizure free post surgery then the EZ has been resected otherwise not. Therefore, we want to build a model to predict the surgical outcome of patients. The first and the most important step in this model building is to design a good set of features that truly encapsulates all the necessary information about EZ. In the following sections, we will be reasoning our way to develop such a set of features.

[0273] 3.4.1 3 dimensional feature vector

[0274] Complete resection of SOZ has been correlated with the patient gaining seizure freedom after surgery. Therefore, the first feature we consider is the extent of SOZ resection. It is captured by the categorical variable, catSOZ, defined as followsI 1. all channels predicted as SOZ bv inter-ictal SOZ classifier is resected ■■ 'f / SOx = <[ 0, odievwise.

[0275] Since our interictal SOZ classifier returns a probability score for each channel being in SOZ, so we can also use the scores to construct a set of features rather than the hard predicted labels as in the catSOZ construction. We construct the second feature by computing what fraction of the total probability score is captured by the resected SOZ nodes. It is captured by the variable, resSOZ ratio, defined as follows

[0276] In equation 3.12, the numerator captures the total probability scores of all the channels that were resected and predicted as SOZ by our interictal SOZ classifier and the denominator captures the total probability scores of all the channels. Similarly, we construct the third feature by computing what fraction of the total probability score is captured by the resected nodes. It is captured by the variable, resection ratio, defined as follows

[0277] In equation 3.13, the numerator captures the total probability scores of all the channels that were resected and the denominator captures the total probability scores of all the channels. Concatenating the 3 constructed features, we get our 3 dimensional feature vector exit SOZjvec resSOZ ratio resection ratio

[0278] 3.4.1.1 Training and Testing corpus for surgical outcome prediction

[0279] We train and validate the performance of the 3 dimensional surgical outcome prediction model on the same multi center patient corpus of 159 patients that we introduced earlier (Table 3.1). Since the patient corpus is unbalanced with ~ 70% patients being seizure free and ~ 30% patients being non-seizure free, so we use 10 fold stratified cross validation to evaluate the performance and generalization power of our model. In 10 fold stratified cross validation, we split the corpus of 159 patients into 10 patients groups of roughly equal size where each group has around 70% seizure free and around 30% non-seizure free patients. Then for each fold, we train the model using 9 patient groups and then test the model on the left out patient group.

[0280] 3.4.1.2 Surgical outcome prediction results using 3 dimensional feature vector

[0281] Recall, from our surgical outcome prediction model setup shown in FIG. 38, that it outputs a probability of the patient being seizure free. Based on the stratified cross validation approach described in the previous section, each patient gets to be in the test set once and hence gets predicted a probability of seizure freedom by the trained model. In FIG. 39, we have plotted the predicted probabilities categorized by the ground truth surgical outcomes. From the plot, it can be observed that the mean predicted probability of seizure freedom is slightly higher for the seizure free group compared to the non-seizure free group.

[0282] We can threshold the probability scores to get the predicted surgical outcome labels and compare them with the ground truth surgical outcome labels. In Table 3.3, we have tabulated the 10 fold performance metric statistics. Since the labels are unbalanced, so Fl score is more relevant than accuracy. The 3 dimensional model is able to predict the surgical outcomes of patients with a mean Fl score of 86% across the 10 folds and also generalizes well as shown by the small variance in performance across folds. Going into the misclassification statistics, we observe that the model makes more mistakes in the minority class, that is the non-seizure free patients, as shown by the 32 false positives.Table 3.3 10 fold performance metric statistics (3 dimensional model): The table presents the aggregated statistics, including the mean and variance, calculated across 10 folds, for four performance metrics: Fl score, Accuracy, Precision, and Recall. The 3 dimensional model demonstrates proficiency in predicting surgical outcomes, achieving a mean Fl score of 86% and a mean accuracy of 79%. Notably, the model exhibits strong generalization capabilities, as evidenced by minimal variance in performance across folds.

[0283] 3.4.1.3 Surgical outcome prediction results using doctor’s SOZ labels

[0284] The gold standard for doctors is to identify the SOZ channels by analyzing all the ictal segments captured in two weeks of hospital monitoring and then surgically resect them with the hope that the patient will become seizure free. We can evaluate the efficacy of the doctors’ process in our context of predicting the post surgical outcomes by defining a categorical variable, docSOZ. which is defined as follows„ I 1, ll channels labeled as SOZ bv doctors are resected docS Z <I 0, otherwise.

[0285] Using the same 10 fold stratified cross validation approach described in Section 3.4. 1. 1, we train and test the 1 dimensional model currently used by the doctors. We observe that the 3 dimensional model highlighted in yellow performs better than the one currently used by the doctors which is highlighted in green, with an accuracy boost of more than 12%, an Fl score boost of more than 7%, and 12 fewer misclassifications.

[0286] 3.4.1.4 Discussion

[0287] We have constructed a 3 dimensional feature vector from the interictal SOZ classifier model scores to predict the surgical outcomes of patients. It performs better than the onecurrently used by doctors and is therefore a better biomarker of EZ than the docSOZ based biomarker used by the doctors currently. We believe that the performance boost is mainly due to two reasons: 1) our interictal SOZ classifier extracts a portion of the unobserved part of EZ which we call PSOZ and 2) the use of the probability scores to construct the features rather than just using the hard predicted labels of the SOZ classifier.

[0288] However, our current 3 dimensional model still makes incorrect predictions on 41 patients. To get a better sense of why our current model make incorrect predictions on these patients, we have plotted the resSOZ ratio as a function of resection ratio for all the 159 patients in FIG. 40. From the plot, it can be observed that there is a subset of seizure free patients (marked with blue circles) which has a low resection ratio and catSOZ = 0 meaning that some of the predicted SOZ were not resected by the doctors. Hence, our current 3 dimensional feature vector (catSOZ, resSOZ ratio, resection ratio) is not able to explain the outcome of these patients because of the low feature values and we suspect that there might be some interesting patterns embedded in the network flows. In order to extract such patterns and augment them into the feature vector, we develop a framework for investigating the network flow dynamics in the next section.

[0289] 3.4.2 5 dimensional feature vector

[0290] In FIG. 40, we found presence of a set of patients who became seizure free even though the entire predicted SOZ was not resected. This evidence led us to investigate the flow dynamics in the constructed interictal synchronization networks. In this framework, we partition the interictal network into 3 non-overlapping regions based on the interictal SOZ classifier predictions and the resection information: predicted SOZ channels that are resected (resected SOZ), predicted SOZ channels that are not resected (preserved SOZ) and predicted Non-SOZ channels that are not resected (preserved Non-SOZ). Based on this partitioning, we compute the flow between these regions by summing the weights of all the edges connecting the nodes in these regions (FIG. 41): fi := Flow between preserved SOZ and preserved Non-SOZ f2:= Flow between resected SOZ and preserved Non-SOZ f? := Flow between preserved SOZ and resected SOZ, where

[0291] Note that such a partitioning might not always exist, when there is no preserved SOZ, in such cases we set fi = f2 = fs = 0.

[0292] In FIG. 42, we have plotted fi, fi, fs as a function of time for a seizure free patient. From the plot, it can be observed that given the resected soz cut, on average there is minimal flow between the preserved SOZ and the preserved Non-SOZ, showing that preserved SOZ cannot trigger large scale synchronization events. Similarly, in FIG. 43, we have plotted fi, f2, fs as a function of time for a non-seizure free patient. In this case, given the resected SOZ cut, there is moderately strong flow between preserved SOZ and preserved Non-SOZ, showing that preserved SOZ channels can trigger large scale synchronization events. Although, we have just plotted the flow evolution for only one seizure free and one non-seizure free patient but such patterns hold for more seizure and non-seizure free patients where such a partitioning exists. Based on these observations, we construct the additional features Flow 1 and Flow 2 which are defined as follows

[0293] In terms of Flow 1 and Flow 2,Seizure free patients := Flow 1 « Flow 2Non-Seizure free patients := Flow 1 ~ Flow 2

[0294] Augmenting the previously constructed 3 dimensional feature vector with the flow features, we get our 5 dimensional feature vector~ [catSOZ, resSOZ ratio, resection ratio, Flow l,Flow 2| '

[0295] We train and validate the performance of the 5 dimensional surgical outcome prediction model using the same 10 fold stratified cross validation approach described in Section 3.4. 1.1.

[0296] 3.4.2. 1 Surgical outcome prediction results using 5 dimensional feature vector

[0297] Recall, from our surgical outcome prediction model setup shown in FIG. 38, that it outputs a probability of the patient being seizure free. Based on the stratified cross validation approach described in Section 3.4.1. 1, each patient gets to be in the test set once and hence gets predicted a probability of seizure freedom by the trained model. In FIG. 44, we have plotted the predicted probabilities categorized by the ground truth surgical outcomes. From the plot, itcan be observed that the mean predicted probability of seizure freedom is higher for the seizure free group compared to the non-seizure free group showing the efficacy of the 5 dimensional feature vector in predicting surgical outcomes.

[0298] We can threshold the probability scores to get the predicted surgical outcome labels and compare them with the ground truth surgical outcome labels. We have tabulated the 10 fold performance metric and misclassification statistics for both the 3 dimensional and 5 dimensional model. The 5 dimensional model which incorporates the network flow dynamics performs better than the 3 dimensional model with a boost in accuracy of 5%, boost in Fl score of 3% and 14 fewer misclassifications. By incorporating the flow dynamics, the model not only corrects the 5 circled Seizure free patients in FIG. 40, but also corrects 9 Non-Seizure free patients.

[0299] 3.4.2.2 Discussion

[0300] We have constructed a 5 dimensional feature vector from the interictal SOZ classifier model scores and the interictal network flow dynamics to predict the surgical outcomes of patients. It performs significantly better than the one currently used by doctors (accuracy boost of 17%, Fl score boost of 10%, 26 fewer misclassifications). Therefore, we have developed an interictal network based method that can predict the post surgical outcome of patients better and hence is a better biomarker of epileptogenic zone than the SOZ based biomarker used by doctors currently. Also, with the interictal network flow dynamics framework we have provided evidence that complete resection of the EZ is not a necessary condition for seizure freedom and sometimes by just resecting the right set of EZ nodes one can disrupt the flow in epileptic network and hence achieve seizure freedom.

[0301] More importantly this framework will allow us to come up with algorithms that can guide the doctors to design minimum resection margins for seizure freedom. For example, in FIG. 45, one such resection margin is shown which will disrupt the flow between the preserved EZ and the preserved Non-EZ resulting in the termination of seizures. The algorithm will have the capacity to generate more resection margin strategies that will disrupt the flow. The thresholds for determining whether a suggested resection plan is good will be determined by retrospective datasets.

[0302] However, our current 5 dimensional model still makes incorrect predictions on 27 patients and the bulk of these mistakes are false positives. To get a better sense of why our current model make incorrect predictions on these patients, we have plotted the resSOZ ratio as a function of resection ratio for all the 159 patients in FIG. 40. Almost all of the 23 falsepositives are located in the far right comer in the dashed box. These patients have a very high resection ratio (> 0.9) and have all the predicted SOZ removed (flow based partitioning doesn’t exist) so the current set of features is not able to explain their outcomes. Therefore, we need to look for some other sources of features that are not necessarily confined to the space of interictal networks. In the next section, we develop such a set of features that combines the interictal network features along with the existing biomarkers of EZ.

[0303] 3.4.3 Mini Seizures: 7 dimensional feature vector

[0304] One of the key observations that led to the accurate identification of EZ from interictal segments was the emergence of hypersynchronous states in the interictal networks. From FIG. 46, we can observe excessive synchronization in the time domain as well during such episodes of hypersynchronization in the interictal period. The observed excessive time domain synchronization in the interictal period looks similar to the ones observed during seizure episodes except it’s less pronounced in terms of amplitude. We hypothesize that these hypersynchronous states are like ’Mini Seizures’ happening in the interictal period with the SOZ / PSOZ nodes driving them and our model is able to mine them for accurate EZ identification.

[0305] Previously, we crafted a weakly supervised deep learning framework aimed at enhancing the specificity of the channel level intenctal biomarker known as HFO by isolating the patho- logical HFO cluster. Our findings indicated that the pathological HFO cluster serves as a superior biomarker of the epileptogenic zone (EZ) compared to the general HFO cluster.

[0306] In Section 3.4.2, we introduced an interictal network property based biomarker of the epileptogenic zone (EZ), which demonstrated superior performance compared to the currently utilized EZ biomarker based on the seizure onset zone (SOZ). To gain insight into the potential relationship between these network level and channel level biomarkers, we generated a plot depicting the time evolution of pathological HFO (spk-HFO) rate of a channel (FIG. 47, left) and time evolution of its eigenvector centrality (FIG. 47, right). The plot suggests a correlation between the eigenvector centrality, representing an interictal network property, and the spk- HFO rate, a channel level interictal biomarker. This observation suggests that channel level pathological events like spk-HFO may be influenced by interictal network effects, implying that we could potentially predict certain statistics of these events using interictal network properties. In the next section, we develop a regression based framework to explore the correlation between spk-HFO and interictal network properties.

[0307] 3.4.3. 1 Regression model

[0308] We built a regression model based framework to capture the correlation between the inter- ictal network properties and spk-HFO rate. Specifically, we predict the smoothed spk- HFO rate of a channel c at time step t using three properties of channel c in the interictal network at time step t: in degree of channel c, out degree of channel c and eigenvector centrality of channel c.

[0309] At a channel level, we quantify the correlation using R2which measures how much of the variations in spk-HFO rate of a channel c can be explained by the interictal network properties.where

[0310] spc(t): true smoothed spk-HFO rate of channel c at time step t

[0311] spcA(t): predicted smoothed spk-HFO rate of channel c at time step t

[0312] E[spc]: Average smoothed spk-HFO rate of channel c

[0313] Rc2: R2score of channel c

[0314] A R2score close to 1 implies that the independent variables (interictal network properties) explain most of the variations in the dependent variable (smoothed spk-HFO rate). An iEEG recording from a patient will have many channels, so to get a patient level feature we summarize the R2 values across channels using mean and the ratio of variance to mean featureteature

[0315] Augmenting the previously constructed 5 dimensional feature vector with the R2 statistics based features, we get our 7 dimensional feature vectorresection ratio j F vec — . Flow 1Flow 2

[0316] We train and validate the performance of the 7 dimensional surgical outcome prediction model using the same 10 fold stratified cross validation approach described in Section 3.4.1.1.

[0317] 3.4.3.2 Surgical outcome prediction results using 7 dimensional feature vector

[0318] Recall, from our surgical outcome prediction model setup shown in FIG. 38, that it outputs a probability of the patient being seizure free. Based on the stratified cross validation approach described in Section 3.4. 1. 1, each patient gets to be in the test set once and hence gets predicted a probability of seizure freedom by the trained model. In FIG. 48, we have plotted the predicted probabilities categorized by the ground truth surgical outcomes. From the plot, it can be observed that the mean predicted probability of seizure freedom is significantly higher for the seizure free group compared to the non-seizure free group showing the efficacy of the 7 dimensional feature vector in predicting surgical outcomes.

[0319] We can threshold the probability scores to get the predicted surgical outcome labels and compare them with the ground truth surgical outcome labels. In Table 3.4, we have tabulated the 10 fold performance metric and misclassification statistics for all the models (1 dimensional, 3 dimensional, 5 dimensional, 7 dimensional). The 7 dimensional model which combines both the interictal network biomarker and channel level pathological event based biomarker performs significantly better than the SOZ based biomarker used by doctors currently with a boost in accuracy of 30%, boost in Fl score of 20% and 50 fewer misclassifications. Also by combining the biomarkers, we have reduced the number of false positives from 23 to only 2.Table 3.4 10 fold performance metric and misclassification statistics compari- son (Idim vs 3dim vs 5dim vs 7dim): The 7 dimensional model significantly outperforms the current 1 dimensional model utilized by doctors, showcasing a substantial increase in accuracy by 30%, a remarkable improvement in Fl score by 20%, and a reduction of 50 misclassifications.

[0320] 3.4.3.3 Discussion

[0321] We have developed a set of features that uses the interictal network properties and exploits the interplay between the network properties and spk-HFO to predict the post surgical outcome of patients better. The significant performance boost in retrospective validation makes the set of features a novel biomarker of epileptogenic zone. Also, the remarkable reduction in invasive monitoring time required to compute the set of features makes this biomarker more appealing as it can considerably cut down the risk and costs associated with hospital stays.

[0322] 5.5 A Biomarker of EZ agnostic surgical resection and SOZ labels

[0323] In the previous section, we observed that adding the features capturing the correlation between interictal network properties and spk-HFO rate significantly boosted the performance of the surgical outcome prediction model. Therefore, there might be some interesting pattern embedded in the correlation feature space that’s driving such good performance. In FIG. 489, we have plottedas a function offor all the 159 patients. Four clusters evolve in this 2 dimensional feature space with two of them predominantly containing seizure free patients and the other two containing non seizure free patients.

[0324] Performing a k-nearest neighbor clustering with 10 fold stratified cross validation on the correlation feature space we get a mean Fl score of 93% (FIG. 50 and Table 3.5) Therefore, there exists an interesting structure in the correlation space that’s enabling the 7 dimensional model to have a near perfect (Table 3.4) surgical outcome prediction power.

[0325] One interesting thing to note is that these two features were constructed by exploiting the correlation between the interictal network properties and the spk-HFO rate and is completely agnostic to surgical resection and SOZ labels. Although, we have unearthed a minimal set of features that can guide surgical resections and help achieve better outcomes in a completely unsupervised manner but there is more information encoded in the four clusters that we observed in FIG. 49. In the next section, we analyze the four clusters and pinpoint some findings that can help the doctors sample the brain regions for optimal surgical resection.Table 3.5 ^-nearest neighbor clustering in the correlation feature space: In the 2-dimensional correlation feature space, we conducted a ^-nearest neighbor clustering with k = l. The clustering algorithm was trained and validated using 10-fold stratified cross- validation, resulting in an impressive mean Fl score of 93%.

[0326] 3.5.1 Analyzing the four clusters

[0327] Clustering the patients in the 2 dimensional correlation space (FIG. 51), we get the following four clusters:

[0328] Cluster 1 : contains non-seizure free patients only and most patients in this cluster have a high mean>|- ' and a high variance in

[0329] Cluster 4: contains mostly seizure free patients and majority of the patients in this cluster have a high mean°-3) and a low variance in

[0330] Cluster 2: contains seizure free patients only and majority of the patients in thiscluster have a low mean < 0-d) and a high variance in ^1

[0331] Cluster 3: contains mostly non-seizure free patients and majority of the patients in this cluster have a low mean<0^) and a low variance in

[0332] 3.5.1. 1 Cluster 1 : high mean high variance cluster of non-seizure free patients

[0333] To gain insights into why patients in Cluster 1 continue to experience seizures despite surgery, we examined the correlation between biomarkers, specifically interictal network properties and spk-HFO rate. We illustrated this correlation by plotting the channel level R2values for a representative patient from Cluster 1 (see FIG. 52). The plot reveals a substantial number of channels with low J?2values indicating poor explanatory power of interictal network properties for the variability in spk-HFO rate within those channels.

[0334] Given that R2quantifies the degree to which interictal network properties can account for the variance in spk-HFO rate, these findings suggest that the constructed interictal network, based on sampled iEEG electrode data, falls short in elucidating the underlying pathological events recorded in certain channels. This discrepancy hints at the presence of an additional epileptic network that remains unexplored by the current electrode placement, potentially responsible for generating these pathological events.

[0335] While the interictal network biomarker proposed in Section 3.4.2, determined by the interictal SOZ classifier prediction scores and network flow dynamics, adequately accounts for the surgical outcomes of patients in cluster 1 — predicted as non-seizure free due to predicted SOZ channels in the preserved region with adequate flow to the rest of the preserved channels — upon closer examination, it appears these patents were unlikely to attain seizure freedom even if all predicted SOZ channels were resected. This is likely attributable to the presence of an unsampled epileptic network.

[0336] 3.5.1.2 Cluster 4: high mean low variance cluster of seizure free patients

[0337] To uncover the reasons behind the seizure free status of patients in Cluster 4 following surgery, we delved into the correlation between biomarkers, specifically examining interictal network properties and spk-HFO rate. Demonstrating this correlation, we plotted channel level R2values for a representative patient from Cluster 4 (refer to FIG. 53). The plot reveals predominantly high R2values across most channels, resulting in a high mean and low variance in R2overall. Given that nearly all spk-HFO occurrences are encompassed by the sampled network, it is highly likely that this precludes the presence of an unsampled epileptic network. Thus, through precise identification and surgical resection of the sampled epileptic network, one can achieve seizure freedom.

[0338] 3.5.1.3 Cluster 2: low mean high variance cluster of seizure free patients

[0339] To understand why patients in Cluster 2 remain seizure-free post-surgery, we closely examined the relationship between two important biomarkers: interictal network properties and spk-HFO rate. This correlation was illustrated by plotting the channel level R2values for a representative patient from Cluster 2 (see FIG. 54). Due to the limited number of detected spk- HFO in patients from these clusters, many channels display zero R2values, contributing to the overall low mean R2.

[0340] However, among the channels where some spk-HFO were detected, they exhibited high R2values, leading to a substantial variance in R2overall. Since almost all spk-HFO occurrences are accounted for by the sampled network, it’s highly probable that this negates the existence of an unsampled epileptic network. Therefore, by accurately identifying and resecting the sampled epileptic network one can achieve seizure freedom.

[0341] 3.5.1.4 Cluster 3: low mean low variance cluster of non-seizure free patients

[0342] To unravel why patients in Cluster 3 continue to experience seizures post-surgery, we examined the correlation between two key biomarkers: interictal network properties and spk- HFO rate. We illustrated this correlation by plotting the channel level R2values for a representative patient from Cluster 3 (see FIG. 55). Given the low number of detected spk- HFO in patients from these clusters, most channels exhibit zero R2values contributing to the overall low mean R2.

[0343] Nevertheless, among the channels where some spk-HFO were detected, their low R2values suggest the presence of an additional epileptic netw ork that remains unexplored by the current electrode placement. This unseen network could potentially be responsible for generating these isolated pathological events, even after surgery.

[0344] The 5 dimensional model (Section 3.4.2), integrating interictal SOZ classifier prediction scores and network flow dynamics, erroneously labeled patients in Cluster 3 as seizure free due to the removal of all predicted SOZ channels. However, a deeper analysis, incorporating the correlation between interictal network properties and spk-HFO rate, in the 7- dimensional model (Section 3.4.3), reveals an unsampled epileptic network. This refinement significantly improves outcome prediction precision, reducing false positives from 23 to 2 (Table 3.4).

[0345] As used in this specification and the claims, the singular forms “a,” “an,” and “the” include plural forms unless the context clearly dictates otherwise.

[0346] As used herein, “about”, “approximately,” “substantially,” and “significantly” will be understood by persons of ordinary skill in the art and will vary to some extent on the context in which they are used. If there are uses of the term which are not clear to persons of ordinary skill in the art given the context in which it is used, “about” and “approximately” will mean up to plus or minus 10% of the particular term and “substantially” and “significantly” will mean more than plus or minus 10% of the particular term.

[0347] As used herein, the terms “include” and “including” have the same meaning as the terms “comprise” and “comprising.” The terms “comprise” and “comprising” should be interpreted as being “open” transitional terms that permit the inclusion of additional components further to those components recited in the claims. The terms “consist” and “consisting of’ should be interpreted as being “closed” transitional terms that do not permit the inclusion of additional components other than the components recited in the claims. The term “consisting essentially of’ should be interpreted to be partially closed and allowing the inclusion only of additional components that do not fundamentally alter the nature of the claimed subject matter.

[0348] The phrase “such as” should be interpreted as “for example, including.” Moreover, the use of any and all exemplary language, including but not limited to “such as”, is intended merely to better illuminate the invention and does not pose a limitation on the scope of the invention unless otherwise claimed.

[0349] Furthermore, in those instances where a convention analogous to “at least one of A, B and C, etc.” is used, in general such a construction is intended in the sense of one having ordinary skill in the art would understand the convention (e.g., “a system having at least one of A, B and C” would include but not be limited to systems that have A alone, B alone, C alone, A and B together, A and C together, B and C together, and / or A, B, and C together.). It will be further understood by those within the art that virtually any disjunctive word and / or phrasepresenting two or more alternative terms, whether in the description or figures, should be understood to contemplate the possibilities of including one of the terms, either of the terms, or both terms. For example, the phrase “A or B” will be understood to include the possibilities of “A” or “B” or “A and B.”

[0350] All language such as “up to,” “at least,” “greater than,” “less than,” and the like, include the number recited and refer to ranges which can subsequently be broken down into ranges and subranges. A range includes each individual member. Thus, for example, a group having 1-3 members refers to groups having 1, 2, or 3 members. Similarly, a group having 6 members refers to groups having 1, 2, 3, 4, or 6 members, and so forth.

[0351] The modal verb “may” refers to the preferred use or selection of one or more options or choices among the several described embodiments or features contained within the same. Where no options or choices are disclosed regarding a particular embodiment or feature contained in the same, the modal verb “may” refers to an affirmative act regarding how to make or use an aspect of a described embodiment or feature contained in the same, or a definitive decision to use a specific skill regarding a described embodiment or feature contained in the same. In this latter context, the modal verb “may” has the same meaning and connotation as the auxiliary verb “can.”

Claims

Claims1. A system for identifying pathological features in medical data comprising: a computing device comprising: a memory storing instructions and multimodal datasets; a processor configured to access the memory to execute the instructions and, thereby, be caused to: receive one or more datasets from one or more modalities associated with brain activity of a population; identify one or more biomarkers for each of the one or more datasets associated with seizure activity; extract one or more features from the one or more biomarkers; identify cognitive spatial biomarkers for each of the one or more datasets; perform a deep learning (DL)-based classification of the one or more biomarkers based on the one or more features; and a user interface configured to output at least one of the one or more biomarkers or the cognitive spatial biomarkers.

2. The system of claim 1, where the processor is further caused to: obtain outcome data based on the one or more biomarkers identified by the processor; and refine at least one of the one or more biomarkers or the cognitive spatial biomarkers.

3. The system of claim 2, wherein the outcome data includes at least one of post-operative seizure freedom or user-identified biomarkers.

4. The system of claim 2, wherein the outcome data includes stimulation data based on the one or more biomarkers identified by the processor.

5. The system of claim 1, wherein the one or more modalities include at least one of invasive brain scan data or non-invasive brain scan data.

6. The system of claim 5, wherein the one or more modalities include at least one of Magnetic Resonance Imaging (MRI), functional MRI (fMRI), Transcranial MagneticStimulation (TMS), Positron Emission Tomography (PET), scalp electroencephalography (EEG), or intracranial EEG (iEEG).

7. The system of claim 1, wherein the one or more biomarkers includes at least one of high frequency oscillations (HFOs) or a network-based biomarker.

8. The system of claim 1, wherein the seizure activity includes at least one of intenctal activity, seizure onset, or ictal activity.

9. The system of claim 1, wherein the one or more biomarkers are classified using either a supervised DL classification, self-supervised DL classification, or an unsupervised DL classification.

10. The system of claim 1, wherein the user interface is further configured to receive parameter settings from the user.

11. A sy stem for identifying treatment target regions of a patient based on pathological features in acquired patient data comprising: a computing device comprising: a memory storing instructions; a processor configured to access the memory to execute the instructions and, thereby, be caused to: receive one or more features extracted from one or more biomarkers associated with seizure activity identified from one or more datasets from one or more modalities associated with brain activity of a population; obtain input data of a patient’s brain comprising one or more modalities to identify the one or more features in the input data; generate a brain map of the patient based on the one or more features identified in the input data; generate a treatment plan for one or more target regions in the patient’s brain based on the brain map; and a user interface configured to receive the input data and output at least one of the one or more biomarkers, the brain map, or the treatment plan.

12. The system of claim 11, wherein the one or more modalities include at least one of invasive brain scan data or non-invasive brain scan data.

13. The system of claim 12, wherein the one or more modalities include at least one of intracranial EEG (iEEG), scalp electroencephalography (EEG), Magnetic Resonance Imaging (MRI), functional MRI (fMRI), Transcranial Magnetic Stimulation (TMS), or Positron Emission Tomography (PET).

14. The sy stem of claim 11, wherein the one or more biomarkers includes at least one of high frequency oscillations (HFOs) or a network-based biomarker.

15. The system of claim 14, wherein the processor further determines a need for additional patient data based on the identified HFO and network-based biomarker.

16. The sy stem of claim 15, wherein the processor further generates an output informing placement of additional EEG electrodes on the patient based on the need for additional patient data.

17. The system of claim 16, wherein the treatment plan includes resection of the one or more target regions.

18. The system of claim 11, wherein the seizure activity includes at least one of interictal activity, seizure onset, or ictal activity.

19. The system of claim 11, wherein the brain map includes at least one of a map of the eloquent cortex or a likelihood heatmap of one or more epileptogenic zones.

20. The system of claim 11, wherein the user interface is further configured to receive parameter settings from the user.

21. A method for identifying pathological biomarkers in medical data comprising: receiving, by a system comprising a processor, one or more datasets from one or more modalities associated with brain activity of a population; identifying, by the processor, one or more biomarkers for each of the one or more datasets associated with seizure activity; extracting, by the processor, one or more features from the one or more biomarkers; identifying, by the processor, cognitive spatial biomarkers for each of the one or more datasets; performing, by the processor, a deep learning (DL)-based classification of the one or more biomarkers based on the one or more features; andoutputting, by a user interface, at least one of the one or more biomarkers or the cognitive spatial biomarkers.

22. The method of claim 21, further comprising: obtaining, by the processor, outcome data based on the one or more biomarkers identified by the processor; and refining, by the processor, at least one of the one or more biomarkers or the cognitive spatial biomarkers.

23. The method of claim 22, wherein the outcome data includes at least one of postoperative seizure freedom or user-identified biomarkers.

24. The method of claim 22, wherein the outcome data includes stimulation data based on the one or more biomarkers identified by the processor.

25. The method of claim 21, wherein the one or more modalities include at least one of invasive brain scan data or non-invasive brain scan data.

26. The method of claim 25, wherein the one or more modalities include at least one of Magnetic Resonance Imaging (MRI), functional MRI (fMRI), Transcranial Magnetic Stimulation (TMS), Positron Emission Tomography (PET), scalp electroencephalography (EEG), or intracranial EEG (iEEG).

27. The method of claim 21, wherein the one or more biomarkers includes at least one of high frequency oscillations (HFOs) or a network-based biomarker.

28. The method of claim 21, wherein the seizure activity includes at least one of interictal activity, seizure onset, or ictal activity.

29. The method of claim 21 , wherein the one or more biomarkers are classified using either a supervised DL classification, self-supervised DL classification, or an unsupervised DL classification.

30. The method of claim 21, further receiving, by the user interface, parameter settings from the user.

31. A method for identifying treatment target regions of a patient based on pathological features in acquired patient data compnsing:receiving, by a system comprising a processor, one or more features extracted from one or more biomarkers associated with seizure activity identified from one or more datasets from one or more modalities associated with brain activity of a population; obtaining, by the processor, input data of a patient’s brain comprising one or more modalities; identifying, by the processor, the one or more features in the input data; generating, by the processor, a brain map of the patient based on the one or more features identified in the input data; generating, by the processor, a treatment plan for one or more target regions in the patient’s brain based on the brain map; receiving input data, by a user interface; and outputting, by the user interface, at least one of the one or more biomarkers, the brain map, or the treatment plan.

32. The method of claim 31, wherein the one or more modalities include at least one of invasive brain scan data or non-invasive brain scan data.

33. The method of claim 32, wherein the one or more modalities include at least one of intracranial EEG (iEEG), scalp electroencephalography (EEG), Magnetic Resonance Imaging (MRI), functional MRI (fMRI), Transcranial Magnetic Stimulation (TMS), or Positron Emission Tomography (PET).

34. The method of claim 31, wherein the one or more biomarkers includes at least one of high frequency oscillations (HFOs) or a network-based biomarker.

35. The method of claim 34, further determining, by the processor, a need for additional patient data based on the identified HFO and network-based biomarker.

36. The method of claim 35, further generating, by the processor, an output informing placement of additional EEG electrodes on the patient based on the need for additional patient data.

37. The method of claim 36, wherein the treatment plan includes resection of the one or more target regions.

38. The method of claim 31, wherein the seizure activity includes at least one of interictal activity, seizure onset, or ictal activity.

39. The method of claim 31, wherein the brain map includes at least one of a map of the eloquent cortex or a likelihood heatmap of one or more epileptogenic zones.

40. The method of claim 31, further receiving, by the user interface, parameter settings from the user.

Images