A COMPUTER IMPLEMENTED METHOD AND COMPUTER PROGRAM PRODUCTS FOR IDENTIFYING TIME-FREQUENCY FEATURES OF PHYSIOLOGICAL EVENTS

Abstract

A method and computer programs for identifying time-frequency features of physiological events are disclosed. A computer system comprises filtering a set of physiological signals within each one of a plurality of time-frequency windows, obtaining a filtered set for each time-frequency window; calculating, for each time-frequency window, a given feature for the filtered set, each one of the signals of the filtered set having a given feature value, providing for each time-frequency window a set of feature values; and calculating, for each time-frequency window a first quantifier defined as a function of said set of features values and/or a second quantifier defined as a function of an empirical distribution of said set of feature values. The first quantifier can be compared with a first threshold and the second quantifier can be compared with a second threshold. The computing system can further select the time-frequency windows satisfying the first threshold and/or the second threshold.

Claims

1. A computer implemented method for identifying time-frequency features of physiological events to extract spectral features of electrophysical seizure onset patterns and to predict an epileptic focus, the method comprising: a) using a computing system to receive: a time period in which a physiological event occurred; a set of physiological signals associated with the physiological event, each one of the set of physiological signals corresponding to a different spatial location of a body part of a living being; a time-frequency region of interest, the time-frequency region being defined by a minimum and a maximum time instant and a minimum and a maximum frequency; the minimum and maximum time instants being comprised within the time period in which the physiological event occurred and the maximum frequency being lower or equal than a sampling rate of the set of physiological signals; and a plurality of time-frequency windows defined within the time-frequency region of interest; b) filtering each one of the set of physiological signals within each of the plurality of time-frequency windows and obtaining as a result a filtered set of physiological signals for each one of the plurality of time-frequency windows; c) using the computing system to calculate, for each one of the plurality of defined time-frequency windows, a given feature for the filtered set of physiological signals, each one of the filtered set of physiological signals having a given feature value, and providing a set of feature values for each time-frequency window and d) using the computing system to calculate, for each time-frequency window, at least one of: a first quantifier defined as a function of the set of features values and comparing the calculated first quantifier with a given first threshold, the computing system further selecting a group of the plurality of time-frequency windows that satisfy the first threshold; and a second quantifier defined as a function of an empirical distribution of the set of feature values and comparing the calculated second quantifier with a given second threshold, the computing system further selecting a group of the plurality of time-frequency windows that satisfy the second threshold.

2. The method of claim 1, wherein the given feature calculated in step c) defines an intrinsic attribute of each one of the filtered set of physiological signals, the intrinsic attribute including at least a power in band (PIB) of each one of the filtered set of physiological signals within each one of the plurality of time-frequency windows or a mean activation (MA), defined as an instantaneous activation of each one of the filtered set of physiological signals averaged within each time-frequency window, the instantaneous activation being a continuous power expressed as a z-score with respect to a common pre-ictal baseline distribution defined by pooling together all signals' power values within a predetermined band.

3. The method of claim 1, wherein the given feature in step c) defines an attribute that indicates how each one of the filtered set of physiological signals is related to other ones of the filtered set of physiological signals, the attribute including at least one of an average Pearson correlation, an average Mutual Information, a betweenness centrality or a node degree of each signal with respect to all other signals within each time-frequency window.

4. The method of claim 1, wherein the first quantifier comprises at least one of the mean, a standard deviation, a maximum, a global activation (GA), a minimum, or a global inactivation (GI), of the set of feature values, the global activation (GA) being defined as a weighted average of a set of positive feature values, where each one of the set of positive feature values is weighted by itself, and the global inactivation (GI) being defined as a weighted average of a set of negative feature values, where each one of the set of negative feature values is weighted by itself.

5. The method of claim 1, wherein the second quantifier comprises at least one of a Renyi entropy, a Fisher information or an activation entropy (AE), of an empirical distribution of the set of feature values; the activation entropy (AE) being defined as a Shannon entropy of the empirical distribution.

6. The method of claim 3, wherein the first quantifier comprises a global activation (GA), a measure defined as a weighted average of a set of positive feature values each one of the set of positive feature values being weighted by itself, and the second quantifier comprising an activation entropy (AE), the activation entropy (AE) being a measure defined as a Shannon entropy of an empirical distribution of the set of feature values.

7. The method of claim 1, further comprising: Comparing the set of feature values with a given third threshold for each one of the group of the plurality of time-frequency windows satisfying the first threshold and/or the second threshold and defining, for each one of the group of the plurality of time-frequency windows satisfying the first threshold and/or second threshold, a subset of the filtered set of physiological signals, collectively referred to as relevant filtered physiological signals; and accumulating all relevant ones of the set of filtered physiological signals across the time-frequency windows satisfying the first threshold and/or second threshold, thus defining a subset of the set of physiological signals, collectively referred to as relevant physiological signals.

8. The method of claim 1, wherein individual ones of the plurality of time-frequency windows overlap with each other.

9. The method of any of claim 1, wherein none of the plurality of time-frequency windows do not overlap with each other.

10. The method of claim 1, wherein individual ones of the plurality of time-frequency windows have an equal width.

11. The method of any of claim 1, wherein individual ones the plurality of time-frequency windows have a different width.

12. The method of claim 1, wherein individual ones of the plurality of time-frequency windows are nested windows with initial bound fixed at the minimum time instant and with increasing final bound.

13. The method of claim 1, wherein the physiological event is an epileptic seizure or a brain response to a delivered stimulus, and the physiological signals being intracranial electroencephalography (iEEG) or scalp electroencephalography (EEG) signals.

14. The method of claim 1, wherein the first and/or the second threshold is the same for all the time-frequency windows.

15. The method of claim 1, wherein each one of the first threshold and the second threshold are determined using a machine learning algorithm applied on the set of physiological signals with prior information on the time-frequency windows of interest.

16. The method of claim 7, wherein the first threshold, the second threshold and/or the third threshold are determined using a machine learning algorithm applied on the set of physiological signals with prior information on the time-frequency windows of interest and the relevant physiological signals.

17. A computer readable medium containing a computer program product for identifying time-frequency features of physiological events to extract spectral features of electrophysical seizure onset patterns and to predict an epileptic focus, the computer program product comprising program instructions that based on: a time period in which a physiological event occurred; a set of physiological signals associated with the physiological event, each one of the set of physiological signals corresponding to a different spatial location of a body part of a living being; a time-frequency region of interest, the time-frequency region of interest being defined by a minimum and a maximum time instant and a minimum and a maximum frequency, the minimum and maximum time instants being comprised within the time period in which the physiological event occurred and the maximum frequency is lower or equal than a sampling rate of the physiological signals; and a plurality of time-frequency windows defined within the time-frequency region of interest: filter each one of the set of physiological signals within each of the plurality of time-frequency windows and obtaining as a result a filtered set of physiological signals for each one of said plurality of time-frequency window; calculate, for each one of the plurality of defined time-frequency window, a given feature for the filtered set of physiological signals, each one of the filtered set of physiological signals having a given feature value, and providing a set of feature values for each time-frequency window; and Calculate, for each time-frequency window, at least one of: a first quantifier defined as a function of the set of features values and comparing the calculated first quantifier with a given first threshold, the program instructions further selecting a group of the plurality of time-frequency windows that satisfy the first threshold; and/or a second quantifier defined as a function of an empirical distribution of the set of feature values and comparing the calculated second quantifier with to a given second threshold, the program instructions further selecting a group of said plurality of time-frequency windows that satisfy the second threshold.

18. The computer readable medium of claim 17, wherein the program instructions further: compare the set of feature values with a given third threshold for each one of the group of said plurality of time-frequency windows that satisfy the first threshold and/or second threshold, thus defining, for each one of the group of the plurality of time-frequency windows satisfying the first threshold and/or second threshold, a subset of the filtered set of physiological signals, collectively referred to as relevant filtered physiological signals; and accumulate all relevant ones of the set of filtered physiological signals across time-frequency windows satisfying the first and/or second threshold, thus defining a subset of the set of physiological signals, collectively referred to as relevant physiological signals.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0030] The previous and other advantages and features will be more fully understood from the following detailed description of embodiments, with reference to the attached figures, which must be considered in an illustrative and non-limiting manner, in which:

[0031] FIG. 1 is a flow chart illustrating a method for identifying time-frequency features of physiological events according to an embodiment of the present invention.

[0032] FIGS. 2A-2C graphically illustrate an embodiment of the method for identifying time-frequency features of physiological events.

[0033] FIGS. 3A-3C is an exploration of spectral activations in different time-frequency windows of interest done around seizure onset.

[0034] FIG. 4 illustrates the processing steps from iEEG signals to SOZ detection.

[0035] FIGS. 5A-5B illustrate the processing steps of an embodiment of the method to identify spectral changes in brain recordings elicited by the administration of a drug.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0036] FIG. 1 illustrates a general overview of a method for identifying time-frequency features of physiological events. The method is implemented/executed by a computer and is applied to a set of physiological signals associated with a physiological event, where each signal of said set corresponds to a different spatial location of a body part of a living being, either human or animal.

[0037] The general purpose of the method is: 1) to identify the characteristic time and frequency scale of spectral changes in said physiological signals associated with (preceding, following and/or co-occurring with) a particular event for which the time boundaries are (approximately) known, and 2) if desired, to determine its spatial localization, i.e., the region where it occurred. In a particular embodiment as will be detailed later, 1) to identify the spectral properties of ictal onset patterns in electrical signals from the brain upon seizure onset (SO), and 2) when seizure onset is spatially confined, to identify the seizure onset zone (SOZ). In other embodiments, as will be detailed later, 1) to identify the spectral features of patterns elicited in brain electrical signals by the administration of a drug or by localized direct electrical stimulation; and 2) to identify brain regions affected by those changes.

[0038] To that end, besides the physiological signals, as input to the method there is also needed a time period in which the physiological event occurred; a time-frequency region of interest defined by a minimum and a maximum time instant comprised within the time period in which the physiological event occurred and a minimum and a maximum frequency, the latter being limited by each signal's sampling rate (in particular, being lower than or equal to the Nyquist frequency of the physiological signals); and a plurality of time-frequency windows defined on the time-frequency region of interest.

[0039] It should be noted that the plurality of time-frequency windows may overlap or not with each other and may have equal or unequal widths. Moreover, in some cases it might be preferable to use nested time-frequency windows to quantify the accumulation of activity, i.e. windows with initial bound fixed at the cited minimum time instant and with increasing final bound.

[0040] Referring back to FIG. 1, at step 1002, the computer filters the set of physiological signals within each of the time-frequency windows, obtaining a filtered set of physiological signals for each time-frequency window as a result. At step 1003, the computer calculates, for each defined time-frequency window, a given feature for the filtered set of physiological signals; where each signal of the filtered set has a given feature value, providing for each time-frequency window a set of feature values. For example, in an embodiment, the given feature is an intrinsic attribute of each signal of the filtered set of physiological signals, including the power in band (PIB) or the mean activation (MA) of each signal within each time-frequency window. Particularly, the MA is defined as the instantaneous activation of each signal averaged within each time-frequency window, where said instantaneous activation is the continuous power (for instance, obtained via the Hilbert transform) expressed as a z-score with respect to a common pre-ictal baseline distribution defined by pooling together all signals' power values within that band. The PIB can be computed using the Hilbert transform method in narrow bands, the wavelet transform, or the Fourier transform. In another embodiment, the given feature defines an attribute that indicates how each signal of the filtered set is related with regard to the other signals of the filtered set. For example, network measures such as node strength or network centrality built upon within-frequency couplings or cross-frequency coupling can be used.

[0041] Then, the computer for each time-frequency window can calculate a first quantifier defined as a function of the set of feature values (step 1004), such as the mean, the standard deviation, the maximum, the global activation (GA), the minimum or the global inactivation (GI). Particularly, the GA and GI target the largest activations and inactivations in the time-frequency windows, respectively, as will be detailed later. Additionally, for each time-frequency window the computer can calculate a second quantifier, defined as a function of an empirical distribution of the set of feature values (step 1007) such as the Fisher information, the Renyi entropy or the activation entropy (AE), that quantifies the structure or confinement of the set of feature values in the time-frequency windows. It should be noted that in other embodiments both first and second quantifiers can be calculated (i.e. both steps 1004 and 1007 can be executed).

[0042] After step 1004, the calculated first quantifier is compared, step 1005, with a first threshold. Finally, at step 1006, the computer selects only those time-frequency windows satisfying the threshold. The same applies for the second quantifier, as shown at steps 1008 and 1009.

[0043] With reference to FIGS. 2A-2C, therein an embodiment of the proposed method is graphically shown. FIG. 2A graphically depicts (1) a set of physiological signals obtained from 5 sensors; and (2) the time period (highlighted box) in which the physiological event occurred, in this case 0-38 seconds. FIG. 2B illustrates (3) the time-frequency region of interest defined by the user (0-30 seconds, 1-150 Hz); and (4) the plurality of time-frequency windows defined on the time-frequency region of interest (nested windows; 4 exemplary windows are shown as highlighted boxes). On the left-side column (or bottom part) of FIG. 2C the different signals filtered within each time-frequency window are illustrated. On the middle column then it is illustrated how the given feature of each one of the filtered signals is calculated. In this example, the given feature is the MA. For each signal, the MA is defined as the instantaneous activation averaged within each time-frequency window, where said instantaneous activation is the continuous power (for instance, obtained via the Hilbert transform) expressed as a z-score with respect to a common pre-ictal baseline distribution defined by pooling together all signals' power values within that band. On the right-side column (or upper part) of FIG. 2C, it is shown how the quantifier for each window is calculated and compared with a threshold. In this example, the quantifier used is the maximum of all feature values within each window and a common threshold of 100 is used. Windows with max(MA)>100 are selected. The window selection reveals that the physiological event defined by the highlighted box in FIG. 2A is characterized by initial HFOs (107-130 Hz) lasting 100 ms (invisible to the naked eye). This pattern of HFOs disappears when considering longer time periods (0-10 s). The window selection also reveals that the event is characterized by a prominent pattern of discharges in the beta band (12-31 Hz) at a longer timescale of appearance (0-20 s). Note that, in this example, time-frequency windows in the beta band (12-31 Hz) and time spans around 0-5 s were not selected by the system (these windows are not shown in the figure), reflecting that patterns in the beta band appear progressively some seconds after the onset of the event (compare with FIG. 2A).

Application to Intracranial EEG for Seizure Focus Detection:

[0044] Following, a particular embodiment of the proposed method for seizure onset pattern identification and epileptic focus (seizure onset zone, SOZ) prediction will be detailed. This is of particular importance in patients with drug-resistant epilepsy, as the accurate localization of the SOZ is crucial to plan a successful surgery. In this embodiment, the physiological event of interest is the epileptic seizure (and in particular the seizure onset) and the physiological signals are EEG recordings acquired from intracranial electrodes. The system is provided with the seizure onset time, and with a time-frequency region of interest around this time, with the plurality of time-frequency windows to be explored.

[0045] To test the method, ten patients with pharmacoresistant epilepsy were selected. These patients underwent stereotactic-EEG presurgical diagnosis at the Epilepsy Unit from Hospital del Mar, Barcelona, Spain. As part of the validation protocol, the EEG epochs containing a total of 67 seizures were selected, annotated and documented by trained epileptologists from this unit and its associated Epilepsy Research Group, within the Neurosciences program at the Hospital del Mar Medical Research Institute (IMIM), Barcelona, Spain. Seizure onset and termination times of each seizure were independently marked by two epileptologists.

[0046] Patient inclusion was based on the following criteria: a) that the seizure focus had been identified by the epileptologists and b) that ictal onset was confined to a reduced number of contacts. After electrode implantation and monitoring, SOZ was identified by neurologists using visual inspection. Surgical resection or radio-frequency thermocoagulation (RFTC) was planned based on individual SEEG evaluations. Intracranial EEG recordings were obtained using 5 to 21 intracerebral multiple contact electrodes that were stereotactically implanted using robotic guidance. Signals were recorded using a standard clinical EEG system with 500 Hz of sampling rate, except for patient 3, where a sampling rate of 250 Hz was used

[0047] For each seizure the computer selected the marked ictal epoch together with 60 seconds of pre-ictal baseline activity. Artifacted channels were identified by visual inspection and removed prior to data analysis. A band-pass filter (FIR, filter band [1,165] Hz) was used to remove slow drifts and aliasing effects. A notch FIR filter at 50 Hz and its harmonic frequencies was also used to remove the alternate current interference.

[0048] Here, the time-frequency region of interest was defined with the following parameters: from 0 to 30 s after seizure onset (SO) and from 1 to 150 Hz. In order to systematically explore different windows for ictal onset pattern identification and for SOZ detection, the following time-frequency grid was defined for exploration. The frequency spectrum was divided into a set of 10 non-overlapping bands. A set of nested time windows obtained by fixing the left bound at SO and varying the right bound were considered. Specifically, the frequency spectrum was splitted into 10 non-overlapping bands using the following cutting points: 1, 4.2, 8, 12, 31, 50, 73, 88, 107, 130 and 150 Hz. In the time-domain, windows with right bound from 100 ms to 30 s after using steps of 100 ms from 100 ms to 5 s and steps of 1 s from 5 to 30 s were considered (i.e., from SO to SO+0.1 s, from SO to SO+0.2 s, . . . , from SO to SO+4.9 s, from SO to SO+5 s, from SO to SO+6 s, . . . , from SO to SO+30 s).

[0049] For each time-frequency window, the mean activation (MA) (Vila-Vidal et al., 2017) of each recorded signal was computed by the system and used as an intrinsic feature of that signal. The MA quantifies the average instantaneous activation of each targeted brain structure for a pre-defined frequency and time windows of interest. The instantaneous activation is the continuous power of each signal (for instance, obtained via the Hilbert transform) in a specific band expressed as a z-score with respect to a single pre-ictal baseline distribution defined by pooling together all signals' power values within that band.

[0050] This quantity is estimated from the EEG signals using a two-step procedure, as described in the next paragraph.

[0051] Let x(t) denote a single-channel EEG signal, and let's assume that the time-frequency window of interest is defined by the frequency range [f.sub.1, f.sub.2] and by time points t.sub.1 and t.sub.2. In practice, x(t) is a discrete time-series denoted by x[k]=x(t), t=kδ, k=1, . . . , K, where δ is the sampling period and K is the total number of samples contained in the recording. Using this notation, time points t.sub.1 and t.sub.2 are indexed as samples k.sub.1 and k.sub.2, respectively. In the first step, a time-varying spectral activation in the frequency range of interest is obtained using the Hilbert transform method. The signal x[k] is band-pass filtered in non-overlapping logarithmically-spaced narrow frequency bands [f, f+Δf] (Δf=0.1 was used) covering the whole frequency range of interest [f.sub.1, f.sub.2]. Signal power in each narrow band is obtained by squaring the signal envelope (modulus of the analytical signal). Summation of the signal power over all narrow frequency bands is performed to obtain the time-dependent power of each region's EEG signal in the desired frequency of interest [f.sub.1,f.sub.2] (PIB). The resulting values are z-scored with respect to a baseline distribution defined by the power values of all contacts in a non-ictal epoch (in this case, from 60 to 20 seconds before ictal onset) to obtain an activation time series A[k] that can be interpreted as a measure of the power change from the pre-ictal state at any given point. Note that A[k] can take both positive and negative values. In the second step, activations are averaged in the time window of interest to obtain the MA for each recording site. For the precise estimation of time-averaged activations, artifact-induced noise should be first removed. Time-stamps with frequency-specific artifacts simultaneously affecting the majority of EEG signals are detected with a sliding-window analysis (200 samples width, 1 sample step), performed independently in the pre-ictal and ictal epochs. Time windows where the product of the mean correlation and the contact-average signal power is two standard deviations larger than the median are considered as artifacts and are discarded in the subsequent analysis. Finally, for a given recording site, the mean activation (MA) in the frequency range [f.sub.1, f.sub.2] and in the time period [t.sub.1, t.sub.2] is obtained by averaging the activation A[k] over the time window of interest (Vila-Vidal et al., 2017):

[00001]$\mathrm{MA}=\frac{1}{k_2-k_1}\underset{k=k_1}{\overset{k_2}{.Math.}}A\left[k\right]$

[0052] Applying this procedure to each recorded signal within a given time-frequency window, one obtains a set of MA values, one for each physiological signal. By repeating the computation within each time-frequency window of the time-frequency grid, the computing system obtains a set of MA values for each window. The MA of a given signal j in the frequency band [f, f+Δf] will be denoted and computed over a time window spanning from the seizure onset until time t with the following notation: MA.sub.j(f,t).

[0053] For optimal focus detection it must be ensured that there is a hierarchical and selective activation of SOZ contacts only. Central to the proposed approach is the definition of two new quantifiers, namely the GA (Global Activation) and the AE (Activation Entropy), that in this particular case are jointly optimized to find time-frequency windows of interest where ictal activity is maximal with respect to a baseline pre-ictal period, while being spatially confined to a few contacts. FIGS. 3A-3C illustrate how these quantifiers are computed and used to assess the amount of information carried in each window of interest. For each time-frequency window, the GA quantifies the magnitude of the most relevant spectral activations with respect to the pre-ictal basal state. It is defined as the weighted average of the set of MA values over all contacts with positive MA, where each contact's contribution is weighted by its MA value, thus ensuring that most active regions have a higher impact on the final value (FIG. 3A):

[00002]$GA\left(f,t\right)=\frac{\underset{j=1}{\overset{N}{.Math.}}{w}_j*M{A}_j\left(f,t\right)}{\underset{j=1}{\overset{N}{.Math.}}{w}_j},$

[0054] with w.sub.j=MA.sub.j(f,t) if MA.sub.j(f,t)>0, and w.sub.j=0, if MA.sub.j(f,t)≤0.

[0055] For each time-frequency window, the AE characterizes the spatial spread of these spectral activations. It is defined as the entropy of an empirical distribution obtained by defining a number of discrete activation ranges on the set of MA values. First, a histogram is computed over the set of MA values using h bins homogeneously spaced between the minimum and maximum MA values (here, h=10 was used). Probability values for each bin (p.sub.i for i=1, . . . , 10) are found as the fraction of contacts lying within the corresponding MA bin (FIG. 3A). The computing system then computes the Shannon's entropy of this empirical distribution using the formula:

[00003]$\mathrm{AE}\left(f,t\right)=\underset{i=1}{\overset{10}{.Math.}}{p}_i*\log p_i.$

[0056] For each time-frequency window in the exploration grid, the computing system extracted the GA and AE quantifiers from the set of MA values. FIGS. 3B and 3C show the GA and AE quantifiers for one exemplary seizure, respectively.

[0057] FIG. 4 summarizes the processing steps to obtain the optimal time-frequency windows for SOZ detection (referred to as seizure onset windows, SOWs) and the SOZ. Intracranial EEG signals in the peri-ictal period were band-pass filtered in pre-defined bands of interest spanning the whole spectrum (FIG. 4A). Then, for each band, MAs were obtained for all possible time windows of interest (FIG. 4B). For each time-frequency window in the exploration grid, the computer extracted the GA and AE from the set of MA values (FIG. 4C).

[0058] Seizure onset window (SOW) detection was achieved by finding time-frequency windows that maximized the GA under the constraint of low AE to ensure that spectral activations were confined only to a few contacts. All pairs (f, t) were considered by the computing system and two threshold conditions, one per variable, were set. Regarding the first variable, a first or lower threshold was set at the 95-th percentile of the GA distribution. For time-frequency-windows to be considered, they were further required to have a GA above 3 to ensure significant global activations with respect to the pre-ictal state. Additionally, a second or upper threshold of 0.5 was set on the AE. All time-frequency windows satisfying both conditions were preselected as candidates to be SOWs. Finally, for each frequency band the first time windows satisfying the condition were only kept. Finally, for each frequency band only the first set of consecutive time windows fulfilling the required criteria was kept. FIG. 4D shows the selected SOWs for the first seizure of patient 1.

[0059] The procedure described before was sequentially applied to the 67 seizures included in the validation. In each seizure, the SOWs were qualitatively found to pinpoint the characteristic frequency and time windows of the seizure onset patterns. As an example, in the first seizure of patient 1, the algorithm selected the following SOWs: 107-130 Hz during the first hundreds of milliseconds and 12-31 Hz between 10 and 20 s. Regions of the posterior and anterior hippocampi were selected as being inside the SOZ. The inspection of the electrophysiological activity around the seizure onset by epileptologists revealed that the output of the method was qualitatively describing the seizure onset patterns. As seen in FIG. 4E the seizure is initiated at a hippocampal level with rapid discharges (˜110 Hz) of very low amplitude in the first hundreds of milliseconds after seizure onset, combined with an inverted depolarizing wave. These discharges are particularly clear in HP1 and are followed by a drastic decrease in frequency evolving to low-voltage fast-activity (LVFA) at 12-31 Hz that becomes visible 5 seconds after seizure onset and that increases in amplitude as the seizure progresses. In this case, activations identified at 12-31 Hz are in fact a combined effect of LVFA activity (˜30 Hz) together with slow rhythmic spikes (RS) (˜15 Hz) of high amplitude particularly observable between 10 and 20 s.

[0060] In each SOW, the relevant iEEG channels were found by applying a third threshold on the set of MA values. In this case, the third threshold was induced by the threshold applied on the AE quantifier. An upper threshold of 0.5 on the AE quantifier implies that, for each time-frequency windows satisfying the threshold, at least 80% of the contacts lie within the same MA bin. For each SOW, the computing system identified the highly populated MA bin (having at least 80% of the contacts), which was used to define the third threshold. Contacts above this bin were considered to be part of the SOZ. This procedure was repeated for all selected SOWs, and SOZ contacts were accumulated, thus obtaining a single SOZ per seizure (FIG. 4F).

[0061] For each patient, the SOZ was defined by accumulating all seizure-specific detected SOZ regions. Then the SOZ defined by epileptologists was used as the ground truth to assess the performance of the proposed method. For each patient, the sensitivity and specificity of the selection were computed. The average sensitivity of the method across patients was 0.94±0.09 (mean±standard deviation), with an average specificity of 0.90±0.09 (mean±standard deviation). In 7 patients all SOZ regions were identified by the method (sensitivity=1). In the remaining 3 false negative cases (SOZ contacts mistakenly marked as non-SOZ) lied at most 1 contact (i.e. 1.5 mm) apart from true positives (regions correctly marked as SOZ).

[0062] The robustness of the method against particular choices of the thresholds was studied by computing a ROC curve as the thresholds varied. Small fluctuations of the GA threshold (95th percentile) and the AE threshold (0.5) were not found to significantly alter the sensitivity and specificity of the results. As the conditions on GA and AE are relaxed more non-SOZ sites are chosen. Despite being outside the SOZ, it was hypothesized that these sites might have a critical role in sustaining and propagating epileptic activity in the early stages of the seizure.

[0063] The amount of SOZ predictability carried in the pre-ictal activity was also assessed. To do so, time windows spanning from seizure onset into the past, i.e., with final bound at seizure onset and initial bound ranging from 30 s to 0.1 s before seizure onset, with the same spacing as before, were also considered. The average sensitivities and specificities of the method across patients in the pre-ictal period were lower than in the ictal period: 0.77±0.32 (mean±standard deviation) and 0.77±0.12, respectively. Although the method has a better performance in the ictal period, high sensitivities and specificities indicate that the pre-ictal period contains sufficient information for SOZ prediction.

[0064] Finally, to relate the results of the proposed method with validated post-operative information, patients that underwent either surgery (n=5) or RFTC (n=4) were sub-selected and the degree of overlap between the predicted SOZ and the treated areas was quantified. Cross-validation with postoperative outcome showed that the fraction of predicted SOZ regions that were treated tended to be higher in patients attaining seizure freedom after a sufficiently long follow-up period that in those exhibiting symptom relapse after the treatment, being this trend higher in surgery than in RFTC. Yet, the proportion of treated SOZ regions lies in the range 0.25-0.85, highlighting the non-trivial relationship between the SOZ and the epileptogenic zone.

[0065] The method for seizure onset pattern identification and epileptic focus prediction was also tested with 9 pharmacoresistant epileptic patients admitted at Hospital Clinic for intracranial video-EEG monitoring between January 2017 and March 2019, including patients with extra-temporal lobe epilepsy. These patients underwent stereotactic-EEG presurgical diagnosis at the EEG lab from the Epilepsy Unit in Hospital Clinic, Barcelona, Spain. Intracranial EEG recordings were obtained using intracerebral multiple contact electrodes that were stereotactically implanted using robotic guidance. Signals were recorded using a standard clinical EEG system with 1024 or 2048 Hz of sampling rate.

[0066] As part of the validation protocol, the EEG epochs containing a total of 44 seizures were selected, annotated and documented by trained epileptologists from this unit. All seizures were processed by the system using the same configuration and procedure as used with patients from Hospital del Mar without any further adaptation or adjustment. Upon revision with clinicians, the system was qualitatively found to pinpoint the ictal onset patterns in a significant proportion of seizures (73%). A preliminary analysis showed that the system correctly identified the SOZ in 6 of the 9 patients when running blindly.

Application to Scalp EEG for Seizure Focus Detection:

[0067] In another embodiment, the method is used to identify ictal onset patterns in scalp EEG and to pinpoint scalp electrodes where those patterns where first registered. In this case, the physiological event of interest is also the onset of an epileptic seizure, but the physiological signals are EEG recordings acquired from scalp electrodes. The present invention was applied for seizure onset pattern identification in scalp EEG recordings. To this end, 2 patients that underwent non-invasive long-term video-EEG monitoring in 2017 in the Epilepsy Unit at Hospital Clinic, Barcelona, Spain, were selected. Video-EEG was performed using a 64-channel Neurolink 64 inbox-1166A amplifier and recorded at 1024 Hz using NeuroWorks software (Natus Medical Inc.). Superficial electrodes were located using the 10/20 International system, using additional electrodes in the frontotemporal regions according the 10/10 system. Video-EEG monitoring was performed in the epilepsy unit for 5 d, and antiepileptic drugs were reduced when necessary to facilitate seizure occurrence.

[0068] EEG epochs containing a total of 4 seizures (2 per patient) were selected, annotated and documented by trained epileptologists from the unit. All seizures were processed by the system using the same configuration and procedure as used with intracranial EEG, except for the AE threshold that was set to 0.7. In all cases, the system identified initial gamma rhythms upon seizure onset that lasted for some seconds before propagation. In particular, the system showed that in some cases HFOs (>130 Hz) that had not been identified in visual inspection by the clinicians were also present in the first milliseconds after seizure onset.

Application to Intracranial EEG for Detection of Brain Responses Evoked by Direct Electrical Stimulation:

[0069] In another embodiment, the system is used to detect spectral patterns elicited by electrical stimulation applied with invasive or non-invasive electrodes and to identify the contacts where these patterns occur. To test this embodiment, intracranial EEG data collected from an epileptic patient that underwent a session of direct electrical stimulation as part of his/her pre-surgical diagnosis at the Epilepsy Unit of Hospital del Mar, Barcelona, was used. In this analyzed patient, seven electrodes were implanted with a total of 80 channels. The sampling rate of the EEG data was 500 Hz. Specifically, the method was applied to analyze the post-stimulation effect that an electrical stimulation delivered on a specific pair of channels in the frontal cingulate area, and which lasted 25 seconds, had on all recorded channels.

[0070] Artifacted channels were identified by visual inspection and removed prior to data analysis. After removal of the stimulation period, a band-pass filter (FIR, filter band [1,165] Hz) was used to remove slow drifts and aliasing effects in the pre-stimulation and post-stimulations periods independently. A notch FIR filter at 50 Hz and its harmonic frequencies was also used to remove the alternate current interference. The system was fed with the EEG signals. In this case, the event of interest was defined as the response to the electrical stimulation during 30 seconds. Here, the same time-frequency region of interest and plurality of time-frequency windows as in the embodiment described for SOZ detection with patients from Hospital del Mar were used. Specifically, the time-frequency region of interest was defined with the following parameters: from 0 to 30 s after stimulation offset (SO) and from 1 to 150 Hz. The frequency spectrum was divided into a set of 10 non-overlapping bands. A set of nested time windows obtained by fixing the left bound at SO and varying the right bound were considered. Specifically, the frequency spectrum was splitted into 10 non-overlapping bands using the following cutting points: 1, 4.2, 8, 12, 31, 50, 73, 88, 107, 130 and 150 Hz. In the time-domain, windows with right bound from 100 ms to 30 s after stimulation offset (SO) using steps of 100 ms from 100 ms to 5 s and steps of 1 s from 5 to 30 s were considered (i.e., from SO to SO+0.1 s, from SO to SO+0.2 s, . . . , from SO to SO+4.9 s, from SO to SO+5 s, from SO to SO+6 s, . . . , from SO to SO+30 s).

[0071] In this case, the feature used to characterize each EEG signal within each time-frequency window was the average power of each signal in a specific band z-scored with respect to the signal's baseline statistics (i.e., data was demeaned by the baseline mean and then normalized by the baseline standard deviation). As a baseline, the system used a 5-seconds artifact-free pre-stimulation period. This quantity was estimated using the following procedure.

[0072] For each signal and frequency band of interest, the system estimated the time-dependent power using the Hilbert transform. Then, resulting values of each contact were z-scored with respect to baseline statistics for each contact and frequency band independently. Then, for each time-frequency window of interest, the system computed the average z-scored power over the time-window of interest.

[0073] For each time-frequency window, the system computed the GA and selected windows with GA>3, searching for increases in power. The system reported general increases in delta, theta, alpha and beta waves (1-30 Hz) that lasted around 30 seconds after stimulation offset. No effect was found in higher frequencies. Additionally, the system computed the AE and selected windows with GA>3 and AE<0.6, searching for spatially confined activations elicited by the stimulation. The system reported an increase in alpha (8-12 Hz) that lasted around 15 seconds and that was visible not only in the frontal cingulate area, but also in distant contacts from the temporal gyri, pinpointing possible propagation paths through fiber tracts.

Application to Intracranial Local Field Potentials for Detection of Brain Responses Evoked by Drug Administration:

[0074] Following, another embodiment of the proposed technology to identify spectral changes in brain recordings elicited by the administration of a drug will be detailed (illustrated in FIGS. 5A-5C). In this embodiment, the physiological event of interest is the administration of the drug and the physiological signals are EEG recordings acquired from scalp or intracranial electrodes. The system is provided with the time period where the drug is infused, with a time-frequency region of interest after this time, with a plurality of time-frequency windows to be explored.

[0075] In particular, the method was tested in intracranial electrophysiological data from a freely-moving mouse during a pharmacological experiment reported in a previous publication (Gener et al., 2019). The experimental paradigm consisted of 30 min baseline period, 30 min following administration of saline, 1 h following administration of the first drug (agonists and antagonists), and 1 h period following administration of the second drug (antagonists following agonists only). During each experiment, intracranial recordings sampled at 30 Khz were obtained from several electrodes implanted in the prefrontal cortex (PFC) and hippocampus (HPC). A complete description of the experimental design and data recordings is provided in (Gener et al., 2019). To test the detector, the PFC and HPC local field potential (LFPs) were reused in one mouse following the administration of 1 mg of 1-(2,5-dimethoxy-4-iodophenyl)-2 amino propane (DOI, partial 5-HT2A/2CR agonist) as a first drug. In total 3 contacts from the PFC and 3 contacts from the HPC were selected. To obtain LFPs, signals were downsampled to 1 kHz, detrended and notch-filtered to remove noise line artifacts (50 and 100 Hz) with custom-written scripts in Python. A band-pass filter (FIR, filter band [1,250] Hz) was used to remove slow drifts and aliasing effects. As part of the preprocessing, signals were also normalized using the z-score.

[0076] In order to detect the time-frequency scale where changes occurred following the drug administration, a set of time-frequency windows were defined. The frequency spectrum was divided into a set of 17 non-overlapping bands using the following cutting points: 1, 4, 8, 12, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140 and 150 Hz. In the time-domain, a set of nested time windows obtained by fixing the left bound at drug administration and varying the right bound were considered. A total of 120 time-windows with right bound from 30 s to 1 h after drug administration (DA) were considered, increasing the window size in steps of 30 s (i.e., from DA to DA+30 s, from DA to DA+60 s, from DA to DA+90 s, . . . , from DA to DA+1 h). In this case, the feature used to characterize each LFP signal within each time-frequency window was the average power in band (PIB) normalized with respect to the baseline PIB distribution of that signal in that particular frequency band. This quantity was estimated using a slicing window and a multi-taper method as described in the next paragraph.

[0077] For each signal, a spectrogram from 1 to 150 Hz was constructed using a slicing window and a multi-taper method. Slicing windows (not to be confused with time windows of interest) with window step and window length equal to 30 s were used. Within each slicing window, power spectral density (PSD) from 1 to 150 Hz was estimated using discrete prolate spheroidal sequences (DPSS) with standardized half bandwidth set to 5 and orders ranging from 1 to 9 (a total of 9 tapers). FIG. 5A shows the spectrograms for two signals from PFC and HPC. Then, for each frequency band defined above, the PIB was obtained by summating the spectrogram over the frequencies contained in that band. This resulted in a time-dependent PIB per signal and frequency band defined above for the whole recording (sampled at 1/30 Hz). For each channel and frequency band, the PIB time-course was z-scored with respect to baseline statistics (i.e., data was demeaned by the baseline mean and then normalized by the baseline standard deviation). As such, the z-scored PIB reflects the power change of a given signal in a given frequency band elicited by the administration of a drug, with positive (resp. negative) values reflecting increases (resp. decreases) in power with respect to the baseline activity of that contact. FIG. 5B shows the z-scored PIB time-courses for two exemplary bands showing spectral changes upon drug administration (4-8 Hz and 90-100 Hz). Finally, for a given time-frequency window defined in the previous paragraph, the feature of the set of LFP signals was obtained by averaging the z-scored PIB values over the time windows of interest. Note that, with a sampling of 1 PIB value every 30 s, the i-th time window of interest (from DA to DA+30i) contains a total of i PIB samples.

[0078] Applying this procedure to each recorded signal within a given time-frequency window, one obtains a set of average z-scored PIB values, one for each physiological signal. By repeating the computation within each time-frequency window of the time-frequency grid, the computing system obtains a set of PIB values for each window. For the sake of clarity, the average z-scored PIB value of a given signal j in the frequency band [f, f+Δf] and computed over a time window spanning from DA until time t is denoted with the following notation: Z.sub.j(f,t).

[0079] For each time-frequency window, two different versions of the method were used. In the first version, the GA (defined above) computed over Z as the first quantifier was used. When computed over Z, the GA is the average of the set of Z values over all contacts with positive Z, where each contact's contribution is weighted by the magnitude of its own Z value, thus ensuring that signals with a most prominent increase in power have a higher impact on the final value:

[00004]$GA\left(f,t\right)=\frac{\underset{j=1}{\overset{N}{.Math.}}{w}_j*Z_j\left(f,t\right)}{\underset{j=1}{\overset{N}{.Math.}}{w}_j},$

[0080] with w.sub.j=Z.sub.j(f,t) if Z.sub.j(f,t)>0, and w.sub.j=0, if Z.sub.j(f,t)≤0.

[0081] Then, time-frequency windows with GA>1.24 (3rd quantile of the GA distribution) were selected. Within each window satisfying this condition, the method further identified contacts responsible for this global increase in power by applying a threshold on the Z values. Contacts with Z>1.24 were selected.

[0082] For the second version, the GI (Global Inactivation) was introduced, a new first quantifier to compute over Z. The GI is defined as the weighted average of the set of Z values over all contacts with negative Z, where each contact's contribution is weighted by the magnitude of its own Z value, thus ensuring that signals with a most prominent decrease in power have a higher impact on the final value:

[00005]$\mathrm{GI}\left(f,t\right)=\frac{\underset{j=1}{\overset{N}{.Math.}}{w}_j*Z_j\left(f,t\right)}{\underset{j=1}{\overset{N}{.Math.}}{w}_j},$

[0083] with w.sub.j=Z.sub.j(f,t) if Z.sub.j(f,t)<0, and w.sub.j=0, if Z.sub.j(f,t)≥0.

[0084] In this case, time-frequency windows with a GI<−1.7 (1st quantile of the GI distribution) were selected. Within each window satisfying this condition, the method further identified contacts responsible for this global decrease in power by applying a threshold on the Z values. Contacts with Z<−1.7 were selected.

[0085] FIG. 5C shows the quantifiers used in the two different versions of the system: the Global Activation (GA) targets the most prominent increases in PIB regardless of the contact where they occur, while the Global Inactivation (GI) targets the most prominent decreases in PIB regardless of the contact where they occur. Drug administration occurs at 0 min.

[0086] Using this procedure, the system reported the following results. Following drug administration, an increase in power in the band 4-8 Hz was observed in some contacts both of PFC and HPC that lasted until the end of the experiment. Additionally, an increase in 40-50 Hz was reported in HPC during the first 20 minutes after DA, that later evolved to a slower rhythm at 12-20 Hz until the end of the experiment.

[0087] Significantly, a power increase in the gamma range (80-110 Hz) was reported in one contact of the PFC, starting 20 minutes after DA and lasting until the end of the recording, in line with the results described in the research article. This power increase in PFC was accompanied by a power decrease in HPC starting 5-10 minutes after DA and lasting until the end of the recording.

[0088] The present invention has been described in particular detail with respect to specific possible embodiments. Those of skill in the art will appreciate that the invention may be practiced in other embodiments. Further, the system and/or functionality of the invention may be implemented via various combinations of software and hardware, as described, or entirely in software elements. Also, particular divisions of functionality between the various components described herein are merely exemplary, and not mandatory or significant. Consequently, functions performed by a single component may, in other embodiments, be performed by multiple components, and functions performed by multiple components may, in other embodiments, be performed by a single component.

[0089] Certain aspects of the present invention include process steps or operations and instructions described herein in an algorithmic and/or algorithmic-like form. It should be noted that the process steps and/or operations and instructions of the present invention can be embodied in software, firmware, and/or hardware, and when embodied in software, can be downloaded to reside on and be operated from different platforms used by real-time network operating systems.

[0090] The scope of the present invention is defined in the following set of claims.

REFERENCES

[0091] Andrzejak R G, David O, Gnatkovsky V, Wendling F, Bartolomei F, Francione S, et al. Localization of epileptogenic zone on pre-surgical intracranial EEG recordings: toward a validation of quantitative signal analysis approaches. Brain Topography 2015; 28(6): 832-837. [0092] Bartolomei F, Chauvel P, Wendling F. Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG. Brain 2008; 131 (7): 1818-1830. [0093] David O, Blauwblomme T, Job A S, Chabardès S, Hoffmann D, Minotti L, et al. Imaging the seizure onset zone with stereo-electroencephalography. Brain 2011; 134(10): 2898-2911. [0094] Geertsema E E, Visser G H, Velis D N, Claus S P, Zijlmans M, Kalitzin S N. Automated seizure onset zone approximation based on nonharmonic high-frequency oscillations in human interictal intracranial eegs. International Journal of Neural Systems 2015, 25(05): 1550015. [0095] Gener T, Tauste Campo A, Alemany-González M, Nebot P, Delgado-Sallent C, Chanovas J, Puig M V. Serotonin 5-HT1A, 5-HT2A and dopamine D2 receptors strongly influence prefronto-hippocampal neural networks in alert mice: Contribution to the actions of risperidone Neuropharmacology 2019; 158: 107743. https://doi.org/10.1016/j.neuropharm.2019.107743 [0096] Gnatkovsky V, Francione S, Cardinale F, Mai R, Tassi L, Lo Russo G, et al. Identification of reproducible ictal patterns based on quantified frequency analysis of intracranial EEG signals. Epilepsia 2011; 52(3): 477-488. [0097] Gnatkovsky V, de Curtis M, Pastori C, Cardinale F, Lo Russo G, Mai R, et al. Biomarkers of epileptogenic zone defined by quantified stereo-EEG analysis. Epilepsia 2014; 55(2): 296-305. [0098] Liu S, Sha Z, Sencer A, Aydoseli A, Bebek N, Abosch A, et al. Exploring the time-frequency content of high frequency oscillations for automated identification of seizure onset zone in epilepsy. Journal of Neural Engineering 2016; 13(2): 026026. [0099] Murphy P M, von Paternos A J, Santaniello S. A novel HFO-based method for unsupervised localization of the seizure onset zone in drug-resistant epilepsy. In: Engineering in Medicine and Biology Society (EMBC), 2017 39th Annual International Conference of the IEEE. IEEE, 2017. P. 1054-1057. [0100] Varatharajah Y, Berry B M, Kalbarczyk Z T, Brinkmann B H, Worrell G A, and Iyer R K. Inter-ictal seizure onset zone localization using unsupervised clustering and bayesian filtering. In: Neural Engineering (NER), 2017 8th International IEEE/EMBS Conference on. IEEE, 2017. P. 533-539. [0101] Vila-Vidal M, Principe A, Ley M, Deco G, Tauste Campo A, Rocamora R. Detection of recurrent activation patterns across focal seizures: Application to seizure onset zone identification. Clinical Neurophysiology 2017; 128(6): 977-985.