- 1Biomedical Engineering Program, Faculty of Engineering, King Abdulaziz University, Jeddah, Saudi Arabia
- 2Neurology Department, Texas Institute for Restorative Neurotechnologies (TIRN), University of Texas Medical School, Houston, TX, United States
- 3Department of Neurological Surgery, University of Alabama at Birmingham, Birmingham, AL, United States
- 4Vivian L. Smith Department of Neurosurgery, McGovern Medical School at University of Texas (UT) Health Houston, Houston, TX, United States
- 5Biomedical Engineering Department, Arizona State University, Tempe, AZ, United States
- 6Department of Translational Neuroscience, Barrow Neurological Institute, Phoenix, AZ, United States
Introduction: The gold standard for identification of the epileptogenic zone (EZ) continues to be the visual inspection of electrographic changes around seizures’ onset by experienced electroencephalography (EEG) readers. Development of an epileptogenic focus localization tool that can delineate the EZ from analysis of interictal (seizure-free) periods is still an open question of great significance for improved diagnosis (e.g., presurgical evaluation) and treatment of epilepsy (e.g., surgical outcome).
Methods: We developed an EZ interictal localization algorithm (EZILA) based on novel analysis of intracranial EEG (iEEG) using a univariate periodogram-type power measure, a straight-forward ranking approach, a robust dimensional reduction method and a clustering technique. Ten patients with temporal and extra temporal lobe epilepsies, and matching the inclusion criteria of having iEEG recordings at the epilepsy monitoring unit (EMU) and being Engel Class I ≥12 months post-surgery, were recruited in this study.
Results: In a nested k-fold cross validation statistical framework, EZILA assigned the highest score to iEEG channels within the EZ in all patients (10/10) during the first hour of the iEEG recordings and up to their first typical clinical seizure in the EMU (i.e., early interictal period). To further validate EZILA’s performance, data from two new (Engel Class I) patients were analyzed in a double-blinded fashion; the EZILA successfully localized iEEG channels within the EZ from interictal iEEG in both patients.
Discussion: Out of the sampled brain regions, iEEG channels in the EZ were most frequently and maximally active in seizure-free (interictal) periods across patients in specific narrow gamma frequency band (∼60–80 Hz), which we have termed focal frequency band (FFB). These findings are consistent with the hypothesis that the EZ may interictally be regulated (controlled) by surrounding inhibitory neurons with resonance characteristics within this narrow gamma band.
Introduction
At least 80 million people (1% of the population) worldwide are afflicted by epilepsy, one-third of whom are refractory (resistant) to antiseizure drugs (Dua et al., 2006). For selected patients with drug-resistant (refractory) epilepsy, surgical resection or ablation of the epileptogenic zone (EZ), conceptualized as the brain region(s) whose obliteration or disconnection is necessary and sufficient to achieve seizure freedom, is an effective therapeutic option (Schmidt, 2009). However, identifying the EZ pre-operatively is challenging as it requires concordance between multiple imaging and electrophysiological investigations (Mégevand and Seeck, 2018). Only if seizure freedom is achieved after surgery, we can conclude that the EZ was within the resected area (Rosenow and Lüders, 2001; Jehi, 2018). Consequently, surgical success is contingent upon accurate localization of the EZ. To this end, the current standard of care includes electrophysiologic source localization from seizures recorded in the Epilepsy Monitoring Unit (EMU) (Shorvon, 2011; Nemtsas et al., 2017). This typically is first attempted via scalp EEG (phase I evaluation). However, in certain cases, EZ localization requires invasive EEG recordings (phase II evaluation) via depth stereo-EEG electrodes or subdural grids (Gelẑiniene et al., 2008). The yield of a singular EMU evaluation is contingent upon the recording of sufficient number of seizures in the EMU and can be inconsistent. Furthermore, prolonged stay in the EMU not only increases cost but also patient’s morbidity (Asano et al., 2005; Kanchanatawan et al., 2014).
On the other hand, a growing body of evidence suggests that the interictal (seizure free) state itself may provide valuable complementary information about the EZ. In previous studies, the EZ is localized by investigating interictal epileptiform discharges (e.g. spikes) as well as non-epileptiform abnormalities like focal slowing (Brodbeck et al., 2011; Englot et al., 2016; Baldini et al., 2020; Conrad et al., 2020). Also, other features of neural activity, ranging from DC shifts, to skew (Mooij et al., 2020), and kurtosis (Akiyama et al., 2012) of signals in pre-defined frequency bands, to high frequency oscillations (HFOs), have been used for localization of the EZ (Modur et al., 2012; Geertsema et al., 2017; Hyde et al., 2019; Thomschewski et al., 2019; Xiang et al., 2020; Stovall et al., 2021). However, it is unclear to what extent the inclusion of interictal epileptiform or other abnormal activity in specific frequency bands can affect the EZ localization’s accuracy. Indeed, in most studies, such feature selection for mapping the EZ is made a priori. Herein, we overcome these parameter selection biases by implementing a big-data analytical approach using multichannel (median of 170 active intracranial EEG [iEEG] channels), long-term (median of 5 days) recordings per patient from phase II evaluation in the EMU, within a broad frequency range (1–500 Hz). This approach employs power ranking, dimension reduction, and clustering techniques.
The discovery of neuronal inductive impedance elucidated a novel mechanism by which neurons are granted frequency preferences (i.e., resonance) (Cole, 1941; Das et al., 2017; Lee et al., 2019). Within these narrow resonant frequencies, neurons generate intrinsic lasting membrane voltage fluctuations and show a heightened response to inputs (Hutcheon and Yarom, 2000; Pike et al., 2000). Increasingly electrophysiological studies have confirmed the existence of similar neuronal resonance frequencies within the EZ, and its pathological role in generating intrinsic frequency-specific epileptiform neural activity (Walton, 1968; Reyes-Garcia et al., 2018). Since seizures (ictal states) are brief paroxysmal events, epileptogenic foci activity in the much longer interictal state is assumed to be constrained from transitioning to the ictal state by nearby inhibition (Prince and Wilder, 1967; Trevelyan et al., 2006; Tsakalis and Iasemidis, 2006; Schevon et al., 2012; Cammarota et al., 2013; Liou et al., 2018). Based on this premise, we hypothesize that the population activity in the EZ during the interictal state should exhibit more intense (time-wise and amplitude-wise) resonance with inhibitory interneurons than those outside the EZ. Thus, an algorithm to localize the EZ during the interictal state could be developed by identifying brain regions that exhibit sustained maximal activity (power) within a narrow band in the frequency domain of inhibition. Our study investigated the existence of such a frequency band across epilepsy patients that also led to identification of brain sites (i.e., iEEG contacts) within the EZ from the interictal period independent of the presence (or absence) of any particular interictal epileptiform activities.
Materials and methods
Subject selection
Patients with medically refractory epilepsy undergoing iEEG monitoring for presurgical evaluation provided informed written consent for analysis of their data. University of Alabama Institutional Review Board approved this study. From a total of 45 patients enrolled in the study between 2016 and 2021 (NeuroNEM database – see Data Availability), 10 patients (six males and four females with a median age of 26.5 ± 15.9 years) met the following inclusion criteria and were included in the training/testing portion of the study (two more patients eventually met the inclusion criteria in 2022 and were thus included in the small, double blinded testing portion of the study): (1) resective surgery was performed to remove the seizure focus and (2) patients were seizure-free for ≥12 months post-surgery (i.e., Engel’s class I). All patients had stereo-EEG electrodes implanted except one patient (PT110) who had subdural grids. Patient selection for invasive EEG monitoring, electrode implantation strategy, duration of recordings, and tapering of antiseizure medications to record seizures were based solely on clinical need and were determined by the independent team of clinicians. Of the 10 patients, eight were diagnosed with mesial temporal lobe epilepsy (TLE) and 2 with extra-TLE (insula and anterior cingulate); the two additional patients also tested had TLE. The 10 patients’ demographics and clinical details are shown in Table 1. Overall, a total of 1,029 h of interictal period from each of 1,400 iEEG channels were analyzed from the 10 patients. A total of 98 seizures were recorded from these patients.
Data acquisition and pre-processing
Long-term continuous iEEG signals were recorded by a Natus Quantum EEG machine (Natus Medical Inc. Pleasant, CA, USA) with a sampling rate of 2,048 Hz. iEEG channels with visually identified artifact for an extended period (longer than 6 h of continuous iEEG) due to electrode or amplifier problems were discarded from subsequent analysis (Only 9 out of 1,409 channels, in 4 out of the 10 patients, were thus contaminated with artifact and were excluded from further analysis). The iEEG data from each patient were then divided into consecutive, non-overlapping T-second epochs. The data in each epoch were pre-processed by detrending and applying a bandpass filter with cutoff frequencies of 1 and 500 Hz, and a notch filter at 60 Hz and its harmonics. We used T of 3 s in the estimation of the univariate periodogram and 10 s for the estimation of the multivariate auto-spectra measures.
Spectral power techniques
In 1897, A. Schuster proposed the periodogram, a non-parametric method of estimating the spectrum in a noisy environment (Schmidt, 1897). He applied the periodogram to find the “hidden periodicities” of the sunspot phenomena, in 1906 (Schmidt, 1906). A parametric method for this application was developed by G. Yule in 1927 (Udnye, 1927) when he fitted the sunspot’s time series with linear regression (i.e., autoregressive – AR model) to detect periodicities in the data. Both of these main methods for estimating spectral power assume white Gaussian noise either superimposed to the data (periodogram method) or driving a recursive harmonic process (autoregressive method). Although, the two methods are computed differently, they both estimated the same intrinsic frequency in the sunspot phenomena (Porat and Marple, 1988). We used both methods of spectral power estimation from the iEEG to investigate the existence of an intrinsic frequency of iEEG sites within the EZ across the two methods. In particular, we estimated the periodogram via discrete Fourier Transform (DFT) per iEEG channel, and the auto-spectra via multivariate autoregressive (MVAR) model fitted to the multi-channel iEEG data.
Univariate spectral power measure (periodogram; S)
The periodogram was estimated per sliding, non-overlapping 3 s epochs per pre-processed multichannel iEEG record. For each epoch, the DFT for the periodogram was estimated by the short-time Fast Fourier Transform (stFFT) per consecutive 250 ms Hamming windows (M) with 50% overlap (hop size = 125 ms). The estimated stFFT matrix of a signal x(t) is given by X(f) = [X1(f),X2(f),…,Xk(f)] such that the mth element of this matrix is:
Xm(f) is the DFT of windowed x(t) data centered around time mR, g(n) is the Hamming window function of length M, and R is the hop size between the successive DFTs. This was done using the stft function in MATLAB R2019b (MathWorks Inc., Natick, MA, USA). To determine the periodogram, the average of the squared magnitude of the successive DFTs per frequency was finally computed per epoch.
Multivariate autoregressive-based spectral power measure (auto-spectra; SKK)
Let y(n) = [y1(n),y2(n),…,yK(n)]T be a K-dimensional vector at time , where fs is the sampling frequency of the data, and with its components being zero mean time series yi(n). Then, for y(n), a K-variate (i.e., K-dimensional) MVAR model of temporal order p (that is, each present value of y(n) depends on p past values of the observed time series yi(n)), is constructed as:
where the order p is determined through criteria from information theory; in our case, using the final prediction error (FPE) method, we estimated p = 25 and kept it at this value across all epochs and patients. The matrices A(τ) contain the model’s (K×K) coefficients (i.e., model parameters) at lag τ(τ = 1,…,p) and were estimated using minimization of the residual noise ε(n) using the Vieira–Morf partial correlation method (Morf et al., 1978). If the model fits the data well, and assuming that each vector component yi(n) is at least a weakly stationary time series within each epoch here of 10 s, the noise (innovation) vector ε(n) = [ε1(n),…,εK(n)]T follows a multivariate standard white noise process having zero mean and covariance matrix . If we denote the (K×K) identity matrix as IK, the MVAR model can be transformed to the frequency domain, as: E(f) = B(f)y(f), where E(f) is the Fourier transform of the residual noise vector and , that is, B(f) essentially results from the Fourier transform of the augmented matrix A of the coefficients of the model (setting A(0) = IK). Then, assuming that ε(n) is the input signal to the model and y(n) the output signal from the model, the transfer function matrix of the model is H(f)=B−1(f). Taking advantage of the spectral factorization theorem (Gevers and Anderson, 1981), the complex power spectral matrix Sij(f) can be calculated directly from H(f) as:
where, the superscript H represents the Hermitian operator.
The pre-processed multichannel iEEG data within consecutive, non-overlapping epochs of 10 s in duration – short enough for the data to be considered statistically stationary and long enough so that we have enough data points for a reliable estimation of the coefficient matrices in the MVAR model – were fitted by MVAR model to compute the complex spectral matrix Sij(f) from (3). Finally, the modulus of the diagonal elements of the matrix (i.e., the auto-spectral power Skk) were computed per channel K and frequency f.
Statistical framework
Reduction techniques
The periodogram (S) and auto-spectra (Skk) methodology yield huge three dimensional matrix of (channel × frequency × time) per patient. Dimensionality reduction techniques were employed in the following order.
Spatial-domain
We applied our statistical ranking approach in order to reduce the spatial dimensionality of the matrix (i.e., ci × frequency × time) by selecting only the channel ci(t,f) that exhibits maximum power at each frequency and time point.
Time-domain
To decrease the dimensionality in the time domain, we then estimated the percentage of time (PoT) the ci channels were identified per frequency within a larger non-overlapping sliding window of 10 min in duration (i.e., each window was composed of 200 epochs of 3 s each in the case of the periodogram or 60 epochs of 10 s each in the case of auto-spectra) then divided by the total number of epochs (200/60). Thus, we end up with the estimation of PoT values for channels that are spatially maxima in power per frequency component (1–500 Hz) every 10 min of the iEEG recording (i.e., ci × frequency × 10 min non-overlapping windows).
Employing the above ranking and time-domain reduction techniques for the periodogram S and auto-spectra Skk power measures, we show in Figure 1 the channel(s) with the highest rate of occurrence (i.e., highest PoT) of maximum power per frequency component (1–500 Hz) throughout the whole EMU recording period (∼103 h) for patient PT103. If the channels selected by this methodology were also clinically assessed focal channels, we denote them by (*) on the plots. In Figure 2, we display the PoT values of only clinically assessed focal channels of the same patient as a heatmap in the time-frequency grid.
Figure 1. Time-frequency plots of the highest percentage of time (PoT) channels that exhibited maximum power within the non-overlapping 10 min window from the two proposed measures [i.e., periodogram (S) – (top) and auto-spectra (Skk) – (bottom)]. Each channel is denoted with a different color. If this methodology selected clinically assessed focal channels for patient PT103, we denote them by (*) in the colorbar. Dashed black lines indicate the timing of the seizures.
Figure 2. Time-frequency heatmap plot of only focal channels from patient PT103 with the highest occurrence rate [highest (PoT)] for the periodogram (left) and auto-spectra (right). Color depth represents the PoT value. Vertical dashed black lines denote the timing of the seizures.
Frequency-domain
In order to reduce the frequency domain, we need to average across an optimal frequency band that renders accurate localization outcomes. Thus, to statistically identify such a band and then validate it for the localization of the EZ, first we introduced the term “focal frequency band” (FFB) as the set of frequency components at which clinically assessed focal channel(s) exhibit the highest PoT. Second, we applied a nested k-fold cross-validation approach to statistically validate the usefulness of FFB in the localization of the EZ, as follows:
1. Apply the Leave-one-patient-out technique to evaluate the success of the localization algorithm, that is, select one patient to leave out for testing and use the remaining nine patients for training of EZILA (e.g., find the optimal / common FFB across all nine patients) and then average across this FFB and (test) the localization algorithm on the 10th patient and see if the iEEG channels thus identified are within that patient’s clinically assessed EZ. Since there is a total of 10 different patients to leave-out, the following steps in training/testing were repeated 10 times, each time for every patient that was left out at this step.
2. In order to increase the statistical power, we use folding in the training set. In particular, select data from a k-fold of the patients in the pool of nine patients (we selected k = 5) to serve as multiple training datasets.
3. Apply the analysis below for each possible subset (fold). Since the order in each subset of patients is not relevant, the number of subsets is given by the combinatorial formula C(9,5) (i.e., = 9!/5!4! = 126 subsets of 5 patients each.
4. Identify the FFB profile (frequency over time), estimated every 10 min (i.e., the PoT matrix), for each of the five patients in a subset.
5. Find each subset’s FFB by computing the intersection across its patients every 10 min.
6. Reduce the time dimension in FFB by creating a histogram of the overlapped FFBs per subset (see Figure 3A).
7. Take the average histogram (AH) across all 126 subsets, then apply K-medoid clustering algorithms to AH, via the kmedoids function in MATLAB R2019b (MathWorks Inc., Natick, MA, USA), to select the range of the FFB (see Figure 3B). Figure 4 displays the FFB derived from the above process for each 10th patient that was left out using the auto-spectra measure of power. It is also noteworthy that the intersection of the thus identified FFBs was almost the same by using either the periodogram or auto-spectra measures of power, that is, (64–76 Hz) and (61–82 Hz), respectively.
8. Average over FFB to reduce the dimension of the PoT matrix to (ci × 10 min non-overlapping windows) for testing the “leave-out” 10th patients.
Figure 3. (A) A histogram of the intersected focal frequency band (FFB) per subset. (B) FFB after averaging histograms (AH) and applying the K-medoid clustering technique as the outliers’ statistics to pinpoint the range of the FFB for one of the “leave-out” 10th patients using the auto-spectra (Skk) measure of power. Red bars represent the selected FFB, and it is (61–83 Hz) for this “leave-out” patient.
Clustering technique and candidate score in testing stage
Data clustering is an unsupervised classification method that aims to create groups (clusters) of objects so that objects in the same cluster are alike and objects in different clusters are significantly distinct (Everitt and Spath, 1985; Arthur and Vassilvitskii, 2007). The number of clusters is determined in advance. Therefore, in the testing stage of the 10th patient, we partitioned the channels in the FFB-averaged PoT matrix into two clusters with the goal of classifying the iEEG channels into either EZ or non-EZ candidate clusters, with the EZ candidate cluster containing channel(s) with high PoT values. The two clusters exhibited a significant gap between them. Note that the EZ candidate cluster could be composed of several channels or just one channel. Then, the candidate score, (range between 0 and 1), of each channel in the EZ candidate group is assigned based on its cumulative presence, within either short or long pre-determined periods, that is, during the first hour of the EMU recording or up to the first typical clinical seizure. Figure 5 provides a flowchart of the aforementioned main steps of EZILA.
Figure 5. Flowchart of the EZ interictal localization algorithm (EZILA). Steps in the rounded corner rectangles are per patient, the step denoted by the oval shape requires information from a training set of patients as described in section “Frequency-domain” in the main text, and the final step in the diagonal rounded rectangle represents the main output of the algorithm per patient.
Results
Testing of “leave-out” 10th patient result
The main outputs of EZILA per patient in the cohort of 10 patients were the five-fold leave-one-out nested process with the clustering and scoring computations. Per the “leave-out” patient, the EZ candidate channels and their scores identified from short and long time periods (i.e., the first hour and up-to-first typical clinical seizure), using the periodogram (S) and the auto-spectra (Skk) measures of power, are shown in Figures 6 and 7, respectively.
Figure 6. The output of the EZ interictal localization algorithm (EZILA) per the (10th) patient via the periodogram (S) measure when applied at two time periods: first hour (blue) and up-to-first seizure (brown). Red tick(s) in the channel label axis denote clinically assessed focal channel(s). EZILA was successful in assigning the highest candidate score to one of the clinically assessed focal channel in all the 10 patients at both time periods. Legends: RAH, right anterior hippocampus; RPH/RPHC, right posterior hippocampus; LAH, left anterior hippocampus; LPH, left posterior hippocampus; RAO/RAOF, right anterior orbitofrontal; RPO/RPOF, right posterior orbitofrontal; LAO/LAOF, left anterior orbitofrontal; LPO/LPOF, left posterior orbitofrontal; LAT, left anterior temporal; RMF, right mesial frontal; RA, right amygdala; RPIN, right posterior insula; RPL, right posterior lesion; SG, left temporal Spencer grid; MG, minigrid; PHC, parahippocampus strip; AH, anterior hippocampus; OF, orbitofrontal; PP, posterior parietal; HC, hippocampus; RMC, right mesial cingulate; RMIN, right mesial insula; RaMC, right anterior mid cingulate; RpPC, right posterior posterior cingulate; LFP, left frontal pole; LpPC, left posterior posterior cingulate; AMG, amygdala; LH, left hippocampus.
Figure 7. The output of the EZ interictal localization algorithm (EZILA) per the (10th) patient via the auto-spectra (Skk) measure when applied at two times periods: first hour (blue) and up-to-first seizure (brown). Red tick(s) in the channel label axis denote clinically assessed focal channel(s). EZILA was successful in assigning the highest candidate score to one of the clinically assessed focal channel in all the 10 patients at the long time period. Legends: see Figure 6.
One important contribution of the EZILA is to identify channels from all recorded iEEG contacts that are likely within the EZ. Each spectral measure contributes toward selecting a few iEEG channels that are probably within EZ. The overall mean and standard deviation of the spatial (channel) reduction across patients was (95.78%; ±2.053) and (96.20%; ±1.593) for S and Skk measures, respectively. This means that EZILA characterized no more than 5% of the available contacts to be within the EZ. This metric of spatial reduction performance, as well as the metric of the percentage of correctly identified focal channels from the reduced channels (i.e., EZ candidate group) taken from the second period (i.e., up-to-first clinical seizure) are shown in Table 2 for both spectral measures.
Retrospectively, to further assess the robustness of the EZILA’s performance over time via the periodogram method, each 10 min window of the EZ candidate cluster was classified into one of three classes based on our prior knowledge from clinical assessment of the actual focal channels per patient: Focal time window if only focal channel(s) compose the EZ candidate cluster, non-Focal time window if only non-focal channel(s) are involved, and Mix time window when the EZ candidate cluster is composed of a mixture of both channel types. The outcome of the above analysis is illustrated for all patients in Figure 8. We note substantial presence of the non-Focal window class in PT105, PT132, and PT143. Interestingly, the former two patients were the only extra-TLE patients (PT105 frontal; PT132 insula) in our cohort and the persistent channel in their non-Focal window class corresponds to the ipsilateral to the focus hippocampal structures. The occurrence of high and frequent power in the gamma range of frequencies in the ipsilateral hippocampus even if the focus is extra-temporal may indicate that the “normal” hippocampus is also involved in the epileptogenic network, which may be more pronounced if hippocampus itself is also pathological, as in the rest of our TLE patients. For PT143, the superficial depth electrode contact (#10) of the focal right anterior hippocampus (RAH) consistently shows up in the non-Focal window class while the focus is more mesial in the right temporal lobe. This constitutes a good example of a more localized identification of the EZ that EZILA could provide us with within the broader area of seizure onset that is presently considered as the EZ.
Figure 8. EZ interictal localization algorithm (EZILA) results overtime via non-overlapping 10 min windows per patient using the periodogram measure. Each time-window is classified as either a focal (red), non-focal (blue), or mix (violet) window based on channel(s) involvement.
Since EZILA is developed as an assistive tool for the EZ localization process in a phase II EMU setting, we consider both the Focal and Mix windows as instants of true positives for localizing at least one site within EZ from the interictal period. In Figure 9 we show the overall percentage of how much assistance can our algorithm provides in the localization process of at least one of the clinically assessed focal channels at any 10 min during phase II EMU recordings per patient.
Figure 9. The EZ interictal localization algorithm (EZILA) robustness in % per patient is measured by the probability of successfully identifying at least one channel within the clinically assessed EZ at any 10 min interictal period of the patient’s (phase II) epilepsy monitoring unit (EMU) recording.
Double-blinded testing of patients
We report the results from two patients (PT1 and PT2), who were only recently recognized as Engel class I, and were not included in previous k-fold training and testing. The FFB used was the intersection of the FFBs across the previous 10 patients (see Figure 4), that is, 64–76 Hz for the periodogram and 61–82 Hz for the MVAR-based measure. The EZILA successfully assigned the highest candidate score to channels within the EZ in both patients The total number of iEEG channels are 102 and 252 for PT1 and PT2, respectively. The performance metric values for these two patients were: 97.06% and 98.41% spatial reduction; and 66.67% and 50% of correctly identified focal channels using the periodogram measure, respectively. In the case of the auto-spectra, the spatial reduction was 96.08% and 98.02; and the percentage of correctly identifying focal channel were 50% and 20% for PT1 and PT2, respectively.
Thus, EZILA correctly assigned the highest score to channels in the EZ from interictal iEEG data analysis in all 12 patients at the (phase II) EMU. It is noteworthy that correctly identifying the EZ in 12 out of 12 patients corresponds to a Jeffrey’s 95% confidence interval of (81.47, 100%) (Brown et al., 2001).
Discussion
Our novel algorithm, EZILA, reliably produces a spatiospectral biomarker of sites in the putative EZ from interictal iEEG data. From hundreds of sampled iEEG contacts, EZILA assigns the highest score to epileptogenic brain regions utilizing their most frequently occurring high power in particular gamma frequency band interictally. This resonance frequency band, common in 10 patients with temporal and extra-temporal epilepsy, ranged from 65 to 75 Hz.
Numerous studies have shown abnormal spontaneous gamma activity during the interictal period. In humans, a scalp EEG study has shown that patients with epilepsy exhibit up to seven times higher gamma power compared to healthy controls (Willoughby et al., 2003). A 2011 study quantified gamma band activity from local field potentials (LFPs) in patients with epilepsy during interictal (1–14 h before seizures) and ictal periods and concluded that there are distinctive bursts of gamma activity in all patients. In addition, a 100-fold increase in gamma power was observed in EZ brain regions (Medvedev et al., 2011).
Local field potential amplitude and frequency depend on the proportional contribution of multiple sources and various properties of the neuronal populations in close proximity to the recording intracranial electrode (Buzsaki et al., 2012). The active membrane property (e.g., intrinsic membrane oscillations) of neurons dictates the generation of the LFP signal (Reimann et al., 2013) and can significantly contribute to its power spectral density (Ness et al., 2016). Furthermore, resonant membrane oscillations must occur synchronously in nearby neurons in order to significantly contribute to the recorded LFP signal, which is a feature that occurs most often in inhibitory interneurons (Buzsaki et al., 2012). A prime example of neurons that exhibit such behavior is GABAergic fast-spiking interneurons (FSIs). They can intrinsically generate and maintain gamma activity of 30–90 Hz (Cardin et al., 2009). Moreover, studies have shown that under a constant tonic input (e.g., epileptiform interictal activity), self-inhibiting populations of interneurons and reciprocally connected pyramidal neuron populations can naturally generate gamma-frequency network oscillations (Whittington et al., 1995, 2000; Börgers and Kopell, 2005).
Based on the inhibitory restraint hypothesis, the EZ is surrounded by a powerful inhibitory drive (i.e., surrounding inhibition), the collapse of which yields an inhibition-excitation imbalance and is believed to produce ictal events (Prince and Wilder, 1967; Treiman, 2001; Tsakalis and Iasemidis, 2006; Tsakalis et al., 2006; Chakravarthy et al., 2007, Chakravarthy et al., 2009; Schevon et al., 2012; Alamoudi, 2022). Surrounding inhibition, as a neural mechanism that creates an inhibitory zone around a central core of activation, has also been reported in the human motor system (Sohn and Hallett, 2004b; Beck and Hallett, 2011; Márquez et al., 2018), movement disorders (Sohn and Hallett, 2004a), cognitive studies (Kiyonaga and Egner, 2016; Shi et al., 2021), and other neurodegenerative disorders (Shin et al., 2007; Belvisi et al., 2021). Our findings are consistent with the inhibitory restraint hypothesis. We observed abnormally sustained gamma resonance-like activity within the EZ during the interictal period, which can represent an emergent property of the pathological firing of inhibitory surrounding interneurons. We hypothesize that because the overall homeostasis of inhibitory-excitatory mechanism in epileptogenic brain regions is defective, “resonating” inhibitory neurons contain the activity of seizure generators, thereby maintaining the seizure-free interictal state. The results from research groups that have considered the interictal EZ-localization problem through directed connectivity analysis from iEEG and magnetoencephalography (MEG) also lead to a containment by a control mechanism in the interictal state (Korzeniewska et al., 2014; Krishnan et al., 2015; Vlachos et al., 2017; Narasimhan et al., 2020). Moreover, Grinenko et al. identified a unique spectro-temporal biomarker of the EZ during the transition from interictal to ictal states irrespectively, of the anatomical location of the focus (Grinenko et al., 2018), which also supports the notion of localization of the EZ based on the dynamical characteristic (sustained) and spectral preference (resonance) of some of its inhibitory constituents independent of its anatomical location.
Reproducing our findings in a bigger cohort of patients would further establish the underlined mechanistic framework. It is noteworthy that electrodes are implanted based upon clinical assessment of EZ for each patient. Therefore, the results, extracted from these electrodes inevitably provide an incomplete picture due to partial coverage of brain activity. Future studies within the presented analytical framework and possibly via noninvasive with more global coverage (e.g., high density scalp EEG, fMRI, MEG) could address this concern.
Conclusion
Analysis of interictal intracranial EEG recorded from patients with focal epilepsy can provide insight into EZ localization. Pathologically increased and sustained activity within a narrow gamma frequency band (65–75 Hz) is consistently observed in epileptogenic regions during the interictal state. By leveraging this property, our novel EZILA correctly identifies channels within the EZ and may thus be used as an assistive tool for early EZ localization during phase II evaluation of patients with epilepsy at the EMU.
Data availability statement
The iEEG datasets analyzed in this article are not readily available because of information that could compromise research participant privacy. Requests to access the de-identified datasets in the NeuroNEM database should be directed to LI, TGVvbmlkYXMuSmFzc2VtaWRpc0BkaWduaXR5aGVhbHRoLm9yZw==; SWFzZW1pZGlzQGFzdS5lZHU=.
Ethics statement
The studies involving human participants were reviewed and approved by the University of Alabama (UAB)’s Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.
Author contributions
OA, LI, and SP contributed to the design of the research. OA contributed to the implementation of the research and the analysis of the results. OA and AI wrote the manuscript. LI and SP reviewed and approved the final version. All authors contributed to the article and approved the submitted version.
Funding
This research study was supported by the US National Science Foundation EPSCoR grant RII-2 FEC OIA1632891 under which the NeuroNEM database was developed (PI: LI).
Acknowledgments
We thank Dr. Ioannis Vlachos, Ms. Farnaz Rezaei, Dr. Emilia Toth, and Dr. Noah Hutson for the numerous fruitful discussions that helped improve this study.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
Akiyama, T., Osada, M., Isowa, M., Go, C. Y., Ochi, A., Elliott, I. M., et al. (2012). High kurtosis of intracranial electroencephalogram as a marker of ictogenicity in pediatric epilepsy surgery. Clin. Neurophysiol. 123, 93–99. doi: 10.1016/j.clinph.2011.05.026
Alamoudi, O. A. (2022). Epileptogenic focus localization by novel analysis of lntracranial EEG in the frequency domain’ dissertation. Ruston, LA: Louisiana Tech University.
Arthur, D., and Vassilvitskii, S. (2007). “K-means++: The advantages of careful seeding”, in Proceedings of the annual ACM-SIAM symposium on discrete algorithms. (Philadelphia, PA: Association for Computing Machinery), 1027–1035.
Asano, E., Pawlak, C., Shah, A., Shah, J., Luat, A. F., Ahn-Ewing, J., et al. (2005). ‘The diagnostic value of initial video-EEG monitoring in children–review of 1000 cases’. Epilepsy Res. 66, 129–135. doi: 10.1016/j.eplepsyres.2005.07.012
Baldini, S., Coito, A., Korff, C. M., Garibotto, V., Ndenghera, M., Spinelli, L., et al. (2020). Localizing non-epileptiform abnormal brain function in children using high density EEG: Electric source imaging of focal slowing. Epilepsy Res. 159:106245. doi: 10.1016/j.eplepsyres.2019.106245
Beck, S., and Hallett, M. (2011). Surround inhibition in the motor system. Exp. Brain Res. 2, 165–172. doi: 10.1007/s00221-011-2610-6
Belvisi, D., Giannì, C., Tartaglia, M., Petsas, N., Baione, V., Crisafulli, S. G., et al. (2021). Abnormal motor surround inhibition associated with cortical and deep grey matter involvement in multiple sclerosis. Clin. Neurophysiol. 132, 1151–1156. doi: 10.1016/j.clinph.2021.01.029
Börgers, C., and Kopell, N. (2005). ‘Effects of noisy drive on rhythms in networks of excitatory and inhibitory neurons’. Neural Comput. 17, 557–608. doi: 10.1162/0899766053019908
Brodbeck, V., Spinelli, L., Lascano, A. M., Wissmeier, M., Vargas, M. I., Vulliemoz, S., et al. (2011). Electroencephalographic source imaging: A prospective study of 152 operated epileptic patients. Brain 134(Pt 10), 2887–2897. doi: 10.1093/brain/awr243
Brown, L. D., Cai, T. T., and das Gupta, A. (2001). Interval estimation for a binomial proportion. Stat. Sci. 16, 101–133. doi: 10.1214/ss/1009213286
Buzsaki, G., Anastassiou, C. A., and Koch, C. (2012). ‘The origin of extracellular fields and currents-EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241
Cammarota, M., Losi, G., Chiavegato, A., Zonta, M., and Carmignoto, G. (2013). Fast spiking interneuron control of seizure propagation in a cortical slice model of focal epilepsy. J. Physiol. 591(Pt 4) 807–822. doi: 10.1113/jphysiol.2012.238154
Cardin, J. A., Carlén, M., Meletis, K., Knoblich, U., Zhang, F., Deisseroth, K., et al. (2009). ‘Driving fast-spiking cells induces gamma rhythm and controls sensory responses’. Nature 459, 663–667. doi: 10.1038/nature08002
Chakravarthy, N., Sabesan, S., Iasemidis, L., and Tsakalis, K. (2007). ‘Controlling synchronization in a neuron-level population model’. Int. J. Neural Syst. 17, 123–138. doi: 10.1142/S0129065707000993
Chakravarthy, N., Tsakalis, K., Sabesan, S., and Iasemidis, L. (2009). ‘Homeostasis of brain dynamics in epilepsy: A feedback control systems perspective of seizures’. Ann. Biomed. Eng. 37, 565–585. doi: 10.1007/s10439-008-9625-6
Cole, K. S. (1941). Rectification and inductance in the squid giant axon. J. Gen. Physiol. 25, 29–51. doi: 10.1085/jgp.25.1.29
Conrad, E. C., Tomlinson, S. B., Wong, J. N., Oechsel, K. F., Shinohara, R. T., Litt, B., et al. (2020). Spatial distribution of interictal spikes fluctuates over time and localizes seizure onset. Brain 143, 554–569. doi: 10.1093/brain/awz386
Das, A., Rathour, R. K., and Narayanan, R. (2017). Strings on a violin: Location dependence of frequency tuning in active dendrites. Front. Cell. Neurosci. 11:72. doi: 10.3389/fncel.2017.00072
Dua, T., Boer, H. M. D., Prilipko, L. L., and Saxena, S. (2006). ‘Epilepsy care in the world: Results of an ILAE/IBE/WHO global campaign against epilepsy survey’. Epilepsia 47, 1225–1231. doi: 10.1111/j.1528-1167.2006.00595.x
Englot, D. J., Nagarajan, S. S, Wang, D. D, Rolston, J. D, Mizuiri, D, and Honma, S. M (2016). The sensitivity and significance of lateralized interictal slow activity on magnetoencephalography in focal epilepsy. Epilepsy Res. 121, 21–28. doi: 10.1016/j.eplepsyres.2016.01.009
Everitt, B., and Spath, H. (1985). Cluster dissection and analysis: Theory, fortran programs and examples. J. R. Stat. Soc. Ser. A 148. 285–286.
Geertsema, E. E., Klooster, M. A. V., Klink, N. E. C. V., Leijten, F. S. S., Rijen, P. C. V., Visser, G. H., et al. (2017). Non-harmonicity in high-frequency components of the intra-operative corticogram to delineate epileptogenic tissue during surgery. Clin. Neurophysiol. 128, 153–164. doi: 10.1016/j.clinph.2016.11.007
Gelẑiniene, G., Endziniene, M., Vaiciene, N., Magistris, M. R., and Seeck, M. (2008). Presurgical evaluation of epilepsy patients. Medicina 44, 585–592. doi: 10.3390/medicina44080076
Gevers, M. R., and Anderson, B. D. (1981). Representations of jointly stationary stochastic feedback processes. Int. J. Control 33, 777–809. doi: 10.1080/00207178108922956
Grinenko, O., Li, J., Mosher, J. C., Wang, I. Z., Bulacio, J. C., Gonzalez-Martinez, J., et al. (2018). A fingerprint of the epileptogenic zone in human epilepsies. Brain 141, 117–131. doi: 10.1093/brain/awx306
Hutcheon, B., and Yarom, Y. (2000). Resonance, oscillation and the intrinsic frequency preferences of neurons. Trends Neurosci. 23, 216–222. doi: 10.1016/S0166-2236(00)01547-2
Hyde, D., Johnson, E., Bolton, J., and Rotenberg, A. (2019). Infraslow eeg analysis refines epileptogenic zone localization. 2019 american clinical neurophysiology society (ACNS) annual meeting. J. Clin. Neurophysiol. 36, 471–524. doi: 10.1097/WNP.0000000000000642
Jehi, L. (2018). The epileptogenic zone: Concept and definition. Epilepsy Curr. 18, 12–16. doi: 10.5698/1535-7597.18.1.12
Kanchanatawan, B., Limothai, C., Srikijvilaikul, T., and Maes, M. (2014). Clinical predictors of 2-year outcome of resective epilepsy surgery in adults with refractory epilepsy: A cohort study. BMJ Open 4:e004852. doi: 10.1136/bmjopen-2014-004852
Kiyonaga, A., and Egner, T. (2016). Center-surround inhibition in working memory. Curr. Biol. 26, 64–68. doi: 10.1016/j.cub.2015.11.013
Korzeniewska, A., Cervenka, M. C., Jouny, C. C., Perilla, J. R., Harezlak, J., Bergey, G. K., et al. (2014). Ictal propagation of high frequency activity is recapitulated in interictal recordings: Effective connectivity of epileptogenic networks recorded with intracranial EEG. Neuroimage 101, 96–113. doi: 10.1016/j.neuroimage.2014.06.078
Krishnan, B., Vlachos, I., Wang, Z. I., Mosher, J., Najm, I., Burgess, R., et al. (2015). Epileptic focus localization based on resting state interictal MEG recordings is feasible irrespective of the presence or absence of spikes. Clin. Neurophysiol. 126, 667–674. doi: 10.1016/j.clinph.2014.07.014
Lee, S. H., Urbano, F. J., and Garcia-Rill, E. (2019). The critical role of intrinsic membrane oscillations. Neurosignals 26, 66–76. doi: 10.1159/000493900
Liou, J. Y., Ma, H., Wenzel, M., Zhao, M., Baird-Daniel, E., Smith, E. H., et al. (2018). Role of inhibitory control in modulating focal seizure spread. Brain 141, 2083–2097. doi: 10.1093/brain/awy116
Márquez, G., Keller, M., Lundbye-Jensen, J., and Taube, W. (2018). Surround inhibition in the primary motor cortex is task-specifically modulated in non-professional musicians but not in healthy controls during real piano playing. Neuroscience 373, 106–112. doi: 10.1016/j.neuroscience.2018.01.017
Medvedev, A. V., Murro, A. M., and Meador, K. J. (2011). Abnormal interictal gamma activity may manifest a seizure onset zone in temporal lobe epilepsy. Int. J. Neural Syst. 21, 103–114. doi: 10.1142/S0129065711002699
Mégevand, P., and Seeck, M. (2018). Electroencephalography, magnetoencephalography and source localization: Their value in epilepsy. Curr. Opin. Neurol. 31, 176–183. doi: 10.1097/WCO.0000000000000545
Modur, P. N., Vitaz, T. W., and Zhang, S. (2012). Seizure localization using broadband EEG: Comparison of conventional frequency activity, high-frequency oscillations, and infraslow activity. J. Clin. Neurophysiol. 29, 309–319. doi: 10.1097/WNP.0b013e318262435d
Mooij, A. H., Frauscher, B., Gotman, J., and Huiskamp, G. J. M. (2020). A skew-based method for identifying intracranial EEG channels with epileptic activity without detecting spikes, ripples, or fast ripples. Clin. Neurophysiol. 131, 183–192. doi: 10.1016/j.clinph.2019.10.025
Morf, M., Vieira, A., Lee, D. T. L., and Kailath, T. (1978). Recursive multichannel maximum entropy spectral estimation. IEEE Trans. Geosci. Electron. 16, 85–94. doi: 10.1109/TGE.1978.294569
Narasimhan, S., Kundassery, K. B., Gupta, K., Johnson, G. W., Wills, K. E., Goodale, S. E., et al. (2020). Seizure-onset regions demonstrate high inward directed connectivity during resting-state: An SEEG study in focal epilepsy. Epilepsia 61, 2534–2544. doi: 10.1111/epi.16686
Nemtsas, P., Birot, G., Pittau, F., Michel, C. M., Schaller, K., Vulliemoz, S., et al. (2017). Source localization of ictal epileptic activity based on high-density scalp EEG data. Epilepsia 58, 1027–1036. doi: 10.1111/epi.13749
Ness, T. V., Remme, M. W. H., and Einevoll, G. T. (2016). Active subthreshold dendritic conductances shape the local field potential. J. Physiol. 594, 3809–3825. doi: 10.1113/JP272022
Pike, F. G., Goddard, R. S., Suckling, J. M., Ganter, P., Kasthuri, N., Paulsen, O., et al. (2000). Distinct frequency preferences of different types of Rat hippocampal neurones in response to oscillatory input currents. J. Physiol. 529, 205–213. doi: 10.1111/j.1469-7793.2000.00205.x
Porat, B., and Marple, S. L. (1988). Digital spectral analysis with applications. Englewood Cliffs, NJ: Prentice-Hall.
Prince, D. A., and Wilder, B. J. (1967). Control mechanisms in cortical epileptogenic foci*: “Surround” inhibition. Arch. Neurol. 16, 194–202. doi: 10.1001/archneur.1967.00470200082007
Reimann, M. W., Anastassiou, C. A., Perin, R., Hill, S. L., Markram, H., and Koch, C. (2013). A biophysically detailed model of neocortical local field potentials predicts the critical role of active membrane currents. Neuron 79, 375–390. doi: 10.1016/j.neuron.2013.05.023
Reyes-Garcia, S. Z., Scorza, C. A., Araújo, N. S., Ortiz-Villatoro, N. N., Jardim, A. P., and Centeno, R. (2018). Different patterns of epileptiform-like activity are generated in the sclerotic hippocampus from patients with drug-resistant temporal lobe epilepsy. Sci. Rep. 8:7116. doi: 10.1038/s41598-018-25378-9
Schevon, C. A., Weiss, S. A., McKhann, G. Jr., Goodman, R. R., Yuste, R., Emerson, R. G., et al. (2012). Evidence of an inhibitory restraint of seizure activity in humans. Nat. Commun. 3:1060. doi: 10.1038/ncomms2056
Schmidt, D. (2009). Drug treatment of epilepsy: Options and limitations. Epilepsy Behav. 15, 56–65. doi: 10.1016/j.yebeh.2009.02.030
Schuster, A. (1897). On lunar and solar periodicities of earthquakes. Proc. R. Soc. Lond. 61, 455–465.
Schuster, A. (1906). On the periodicities of sunspots. Philos. Trans. R. Soc. Lond. Ser. A 206, 69–100.
Shin, H. W., Kang, S. Y., and Sohn, Y. H. (2007). Disturbed surround inhibition in preclinical parkinsonism. Clin. Neurophysiol. 118, 2176–2179. doi: 10.1016/j.clinph.2007.06.058
Shi, R., Gao, H., and Zhang, Q. (2021). The extent of center-surround inhibition for colored items in working memory. Mem. Cogn. 49, 733–746. doi: 10.3758/s13421-020-01116-3
Shorvon, S. (2011). Handbook of epilepsy treatment, 3rd Edn. New York, NY: Wiley-Blackwell. doi: 10.1002/9781444325201
Sohn, Y. H., and Hallett, M. (2004b). Surround inhibition in human motor system’. Exp. Brain Res. 158, 397–404. doi: 10.1007/s00221-004-1909-y
Sohn, Y. H., and Hallett, M. (2004a). Disturbed surround inhibition in focal hand dystonia. Ann. Neurol. 56, 595–599. doi: 10.1002/ana.20270
Stovall, T., Hunt, B., Glynn, S., Stacey, W. C., and Gliske, S. V. (2021). Interictal high frequency background activity as a biomarker of epileptogenic tissue. Brain Commun. 3:fcab188. doi: 10.1093/braincomms/fcab188
Thomschewski, A., Hincapié, A. S., and Frauscher, B. (2019). Localization of the epileptogenic zone using high frequency oscillations. Front Neurol. 10:94. doi: 10.3389/fneur.2019.00094
Treiman, D. M. (2001). GABAergic mechanisms in epilepsy. Epilepsia 42(Suppl. 3) 8–12. doi: 10.1046/j.1528-1157.2001.042Suppl.3008.x
Trevelyan, A. J., Sussillo, D., Watson, B. O., and Yuste, R. (2006). Modular propagation of epileptiform activity: Evidence for an inhibitory veto in neocortex. J. Neurosci. 26, 12447–12455. doi: 10.1523/JNEUROSCI.2787-06.2006
Tsakalis, K., Chakravarthy, N., and Sabesan, S. (2006). A feedback control systems view of epileptic seizures. Cybern. Syst. Anal. 42, 483–495. doi: 10.1007/s10559-006-0087-2
Tsakalis, K., and Iasemidis, L. (2006). Control aspects of a theoretical model for epileptic seizures. Int. J. Bifurc. Chaos 16, 2013–2027. doi: 10.1142/S0218127406015866
Udnye, Y. G. (1927). On a method of investigating periodicities in disturbed series, with special reference to wolfer’s sunspot numbers. Philos. Trans. R. Soc. Lond. Ser. A Containing Pap. Math. Phys. Character 226, 267–298.
Vlachos, I., Krishnan, B., Treiman, D. M., Tsakalis, K., Kugiumtzis, D., Iasemidis, L. D., et al. (2017). The concept of effective inflow: Application to interictal localization of the epileptogenic focus from IEEG. IEEE Trans. Biomed. Eng. 64, 2241–2252. doi: 10.1109/TBME.2016.2633200
Walton, J. N. (1968). Epilepsy (contemporary neurology series, no. 2): By RP Schmidt and BJ Wilder, viii+ 220 pages, 49 illustrations, 9 tables, Blackwell. Oxford: Elsevier.
Whittington, M. A., Traub, R. D., Kopell, N., Ermentrout, B., and Buhl, E. H. (2000). Inhibition-based rhythms: Experimental and mathematical observations on network dynamics. Int. J. Psychophysiol. 38, 315–336. doi: 10.1016/S0167-8760(00)00173-2
Whittington, M. A., Traub, R. D., and Jefferys, J. G. R. (1995). Synchronized oscillations in interneuron networks driven by metabotropic glutamate receptor activation. Nature 373, 612–615. doi: 10.1038/373612a0
Willoughby, J. O., Fitzgibbon, S. P., Pope, K. J., Mackenzie, L., Medvedev, A. V., Clark, C. R., et al. (2003). Persistent abnormality detected in the non-ictal electroencephalogram in primary generalised epilepsy. J. Neurol. Neurosurg. Psychiatry. 74, 51–55. doi: 10.1136/jnnp.74.1.51
Keywords: epilepsy, epileptogenic focus, interictal localization, seizure onset zone (SOZ), resonance, surround inhibition, spectral analysis, long-term sEEG
Citation: Alamoudi OA, Ilyas A, Pati S and Iasemidis L (2022) Interictal localization of the epileptogenic zone: Utilizing the observed resonance behavior in the spectral band of surrounding inhibition. Front. Neurosci. 16:993678. doi: 10.3389/fnins.2022.993678
Received: 14 July 2022; Accepted: 16 November 2022;
Published: 12 December 2022.
Edited by:
Paul Carney, University of Missouri System, United StatesReviewed by:
Petros Nemtsas, Aristotle University of Thessaloniki, GreeceAndrei V. Medvedev, Georgetown University, United States
Copyright © 2022 Alamoudi, Ilyas, Pati and Iasemidis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Omar A. Alamoudi, b2FtYWxhbW91ZGlAa2F1LmVkdS5zYQ==
†These authors share senior authorship