Resting-state background features demonstrate multidien cycles in long-term EEG device recordings ================================================================================================= * William K.S. Ojemann * Brittany H. Scheid * Sofia Mouchtaris * Alfredo Lucas * Joshua J. LaRocque * Carlos Aguila * Arian Ashourvan * Lorenzo Caciagli * Kathryn A. Davis * Erin C. Conrad * Brian Litt ## Abstract **Background** Longitudinal EEG recorded by implanted devices is critical for understanding and managing epilepsy. Recent research reports patient-specific, multi-day cycles in device-detected epileptiform events that coincide with increased likelihood of clinical seizures. Understanding these cycles could elucidate mechanisms generating seizures and advance drug and neurostimulation therapies **Objective/Hypothesis** We hypothesize that seizure-correlated cycles are present in background neural activity, independent of interictal epileptiform spikes, and that neurostimulation may disrupt these cycles. **Methods** We analyzed regularly-recorded seizure-free data epochs from 20 patients implanted with a responsive neurostimulation (RNS) device for at least 1.5 years, to explore the relationship between cycles in device-detected interictal epileptiform activity (dIEA), clinician-validated interictal spikes, background EEG features, and neurostimulation. **Results** Background EEG features tracked the cycle phase of dIEA in all patients (AUC: 0.63 [0.56 - 0.67]) with a greater effect size compared to clinically annotated spike rate alone (AUC: 0.55 [0.53-0.61], p < 0.01). After accounting for circadian variation and spike rate, we observed significant population trends in elevated theta and beta band power and theta and alpha connectivity features at the cycle peaks (sign test, p < 0.05). In the period directly after stimulation we observe a decreased association between cycle phase and EEG features compared to background recordings (AUC: 0.58 [0.55-0.64]). **Conclusions** Our findings suggest that seizure-correlated dIEA cycles are not solely due to epileptiform discharges but are associated with background measures of brain state; and that neurostimulation may disrupt these cycles. These results may help elucidate mechanisms underlying seizure generation, provide new biomarkers for seizure risk, and facilitate monitoring, treating, and managing epilepsy with implantable devices. * Background EEG features track multidien cycles in RNS dIEA * dIEA linked EEG features are patient-specific and may differ with cortical structure * Responsive neurostimulation suppresses EEG-dIEA coupling Keywords * neurostimulation * epilepsy * RNS * biomarkers * interictal * spikes ## Introduction Temporal cycles in seizures have been documented extensively in animals and humans for many years [1]–[4]. Recently, recordings from chronic implantable EEG recording devices have yielded new insights into seizure generation in patients with epilepsy [5], [6]. In particular, Baud et al. found patient-specific, periodic multi-day (“multidien”) fluctuations in device-detected interictal epileptiform activity (dIEA) over the course of weeks to months in long-term recordings from the implantable responsive neurostimulation (RNS) device [5]. More recently, similar cycles have been described in long-term EEG recordings from other devices in humans [6]–[8], dogs [9], and rats [10], and in longitudinal measures of heart rate [11]. Device-detected IEA counts are based on combinations of detection parameters that are hand-tuned by neurologists to be sensitive to acute aberrant brain activity. In practice, these detection parameters vary over time, and it is unclear how often they track epileptiform discharges or some other phenomena. Despite the heterogeneous and subjective nature of these detection settings, there is evidence that tracking patient-specific dIEA can be useful as a predictive biomarker of seizure risk [5], [12]. While numerous theories exist regarding the origin of multidien cycles, and early evidence indicates that cycles can be perturbed with neuromodulation [7], little is known about how their biological underpinnings relate to neural function and therapy [4]. Understanding these cycles and their relationship to other measures of brain state would allow us to discover novel ways to titrate treatment and support patient-specific, phase-locked therapeutic paradigms. Prior efforts to investigate multidien fluctuations in long-term neural recordings are based upon measures of acute epileptiform activity (i.e. spikes, dIEA) [5], [8], [13], however it is unclear whether these cycles result from epileptiform activity or some other phenomena. The relationship between dIEA cycles and background EEG measures remains unexplored. In this study we investigate the relationship between interictal background EEG features and dIEA cycle phase. We test the hypothesis that the dIEA cycle that best aligns with periodic seizure occurrence is represented in background neural activity as well as in clinical spike rate, and that neurostimulation suppresses this relationship. To test these hypotheses, we measure associations between dIEA cycle phase and features calculated from interictal recordings: spike rate, band power, and connectivity both at baseline and in a window following stimulation. We identify novel associations between seizure-associated dIEA cycles and measures of background brain state, suggesting that the observed cycles are not solely based on acute discharges and creating a platform for generating novel hypotheses that further our understanding and treatment of epilepsy. ## Materials and Methods ### RNS system recordings The RNS System is a neuromodulation therapy that uses an implantable device with two bipolar channels each in two electrode leads (4 channels total). EEG activity in two selected detection channels is continuously monitored, and the device responds with electrical stimulation pulses when abnormal activity is detected [14]. A count of hourly detections is recorded by the device and has previously been referred to as a measure of interictal epileptiform activity (IEA) [5]. Here we denote device-detected interictal epileptiform activity as “dIEA”, to distinguish it from clinically-verified epileptiform spikes that we extract from EEG recordings (in methods below). In addition to recording counts of hourly dIEA detections, the RNS device records 90-second “Scheduled Event” (SE) EEG recordings every 12 hours to capture baseline activity (**Figure 1D**), and 90-second “Long Episode” (LE) recordings, which capture likely seizures. Due to device memory limitations, SE clips are often set as low priority and can be overwritten by other stored EEG clips. Because the detection criteria are met so frequently, SEs often capture stimulations as well as baseline activity (**Figure 1F**). In this study, we used hourly dIEA detections, LE capture timepoints, and EEG signals recorded in SE events both with and without stimulation. ![Figure 1](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/07/07/2023.07.05.23291521/F1.medium.gif) [Figure 1](http://medrxiv.org/content/early/2023/07/07/2023.07.05.23291521/F1) Figure 1 EEG clip selection and dIEA phase labeling. We illustrate our methods in panels A-E using data from the same representative patient (HUP096). **A)** Hourly counts of detected interictal epileptiform activity (dIEA) after z-scoring between patient visits (vertical red lines). A 100-h moving average was applied for visualization. **B)** The periodogram obtained by averaging the spectrogram of the dIEA signal over time. The peak associated with the greatest LE PLV is indicated by the purple marker. **Inset:** Histogram showing mean multidien phase and axial mean of LE occurrences. LEs in this patient were phase-locked to the rising phase of the multidien cycle. **C)** Signal component after wavelet decomposition of dIEA signal with period corresponding to the marked periodogram peak in **B**. Sawtooth lines indicate instantaneous phase, purple markers indicate LE recordings, filled dots indicate sampling of SE recordings in peak (P), falling (F), trough (T), and rising (R) phase bins. **D)** Example SE recording without stimulation sampled at the cycle trough. Each channel represents a bipolar montage between pairs of adjacent contacts on the same electrode. **E)** Distribution of recorded stimulation-free SEs across each phase bin. **F)** Example SE recording with shaded 10-s analysis windows following simulations marked in blue. **G)** Distribution of post-stimulation analysis windows across each phase bin. LE-long episode, PLV-phase-locking value, SE-scheduled event, R-rising, P-peak, F-falling, T-trough ### Patient population We retrospectively analyzed longitudinal RNS EEG and dIEA recordings from twenty out of 28 patients with drug-resistant epilepsy who were implanted with the RNS System (NeuroPace, Inc., Mountain View, CA) at the Hospital at the University of Pennsylvania. We included patients if they had at least one year of dIEA recordings (recordings could be segmented as long as segments were at least three months in duration). Patients were required to have at least 200 recorded SE clips without stimulation, or 200 post-stimulation analysis windows, so that enough samples were present for our multivariate analysis. We also excluded patients without a significant multidien cycle in their dIEA counts. Eight patients did not meet our inclusion criteria (full breakdown is detailed in the **Supplemental Methods**). The remaining 20 patients included in this study were implanted with the RNS device between August 2015 and January 2022. Patient demographics are summarized in **Table 1**. All patients considered for this study gave written informed consent in accordance with the Institutional Review Board of the University of Pennsylvania. View this table: [Table 1.](http://medrxiv.org/content/early/2023/07/07/2023.07.05.23291521/T1) Table 1. Patient Demographics ### Extracting multidien cycles from dIEA counts We followed a similar methodology as outlined in Baud et al. to extract multidien cycles from dIEA activity (**Figure 1A**) [5], [15]. In brief, for each patient we applied a wavelet transform to the dIEA signal and identified the dominant modes in the average frequency spectrum (peaks in the periodogram) in the range of 3-60 days. We performed a wavelet decomposition at each peak periodicity and for each patient, we used the wavelet component that was most phase-locked to seizure occurrence for our analysis (**Figure 1B, C**). Finally, we calculated phase angle by applying the Hilbert transform to each patient’s extracted multidien cycle (**Figure 1C**). See **Supplemental Methods** for more details. ### Extracting electrophysiology from RNS EEG recordings #### Selecting device EEG recordings We analyzed only SE device recordings because we were interested in understanding the relationship between *baseline* brain state and dIEA cycle phase. We divided our analysis of SEs into two separate datasets: (1) SEs without detections and subsequent simulations (N=19 patients), and (2) SEs containing simulations followed by a 12-second window of stimulation-free EEG (N = 17 patients) (**Figure 1D, F**). SEs are non-uniformly distributed in time (**Figure S13**) precluding us from applying circular analyses and wavelet decompositions to directly compare EEG recordings with dIEA cycles. We assigned a phase angle to each SE based on its phase location within the extracted dominant multidien cycle. The SE clips were categorized into four phase bins: peak (-π/4 to π/4), falling (π/4 to 3π/4), trough (3π/4 to -3π/4), and rising (-3π/4 to -π/4). #### Extracting spike rate We detected spikes in all SE clips using a modified version of a previously validated spike detector [16] (**Figure 2A**) (see **Supplemental Methods** for detector parameters). One board-certified neurologist (JL) independently validated a random selection of 50 spikes per patient. The median positive predictive value (PPV) across patients was 70% (IQR: [47.5%, 78.5%]). ![Figure 2](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/07/07/2023.07.05.23291521/F2.medium.gif) [Figure 2](http://medrxiv.org/content/early/2023/07/07/2023.07.05.23291521/F2) Figure 2 Spike rate association with dIEA counts and multidien cycles. A) Illustration of spike detection limits for an example spike detection. B) Spike rate in scheduled events and dIEA counts for sample subject HUP096. For visualization purposes, dIEA counts were smoothed using a moving average with window size equal to the median schedule event sampling rate (∼12 hours) and the averaged dIEA signal was resampled to correspond with the recorded scheduled events. C) Polar plot for sample subject HUP096 of average spike rate (spikes/90s clip) and standard error across phases of their multidien cycle. D) Box plots show distribution of spike rate in rising/falling (green) and peak/trough (orange). E) Cohen’s D effect sizes of spike rate association in PvT (green) or RvF (orange). Star represent a significant difference in phase groups (signtest). Population distribution of effect sizes i represented in the box plot at right, with the distribution of PvT effect sizes being significantly greater than zero (signtest). #### Extracting representative brain-state features We calculated band power and connectivity in four frequency bands - theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz), and gamma (33-100 Hz) (**Table 2**). In SEs without stimulation, we calculated features across the entire 90s clip, in SEs with stimulation, features were calculated within the 10-second window beginning two seconds after the stimulation occurred (**Figure 1F)** [17]. We calculated the maximum feature value between detection channels (band power) and between leads (connectivity) as a summary statistic for each feature in each window. We also calculated a binary time-of-day feature, delimited at 8am/8pm, to account for circadian variation in the features when analyzing their model coefficients. To address changes in detection parameters and baseline drifts in the features, we z-scored each feature between clinical neurologist visits. See the **Supplementary Methods** for details on feature calculations. View this table: [Table 2:](http://medrxiv.org/content/early/2023/07/07/2023.07.05.23291521/T2) Table 2: Feature calculations. Two types of features were extracted from the EEG recordings measuring (1) neural activity at each detection channel and (2) functional connectivity within leads. Each feature was calculated in four canonical frequency bands encapsulating the entire frequency spectrum captured by the RNS device. For each bandpower feature we took the maximum value between the two detection channels (purple vs green). For each connectivity feature we took the maximum value between the two leads (purple vs green). We also included time of day of the recording to account for potential circadian variations when interpreting our features. ### Statistical analysis #### Multidien cycle significance and phase distribution We tested the significance of the dIEA periodogram peak associated with the highest LE PLV value using the 99% confidence interval of an autoregression-based red noise null model for periodic signals, a method previously applied to find multidien cycles in humans [18], [19]. When analyzing the phase entrainment of LEs to cycle phase, we used an omnibus test for non-uniformity to assess the significance of LE phase locking [5], [20]. #### Univariate feature analysis For each univariate feature extracted from the SEs, we used a one-way ANOVA to test for significant differences in feature values across the four phase bins. We further split the feature comparison into two groups to increase interpretability: peak vs. trough (PvT) phase bins, and rising vs. falling (RvF) phase bins. We calculated the effect size (Cohen’s *d*) between the peak and trough groups, and the rising and falling groups. We report effect size as Cohen’s *d* to better characterize both the directionality and sample-size-independent magnitude of the relationships[21]. We tested for group-level trends in effect-size and directionality across the patient population using a two-tailed sign test against a zero median for each feature (**Figure S5**). #### Multivariate machine learning model We developed a multivariate classification model to understand how band power and connectivity features relate to dIEA cycle phase. For each patient, we built two linear support vector machines (SVMs) which used the EEG features to classify PvT and RvF phases, respectively. We transformed the features using principal component analysis (PCA) to account for multicollinearity. For each patient-specific PvT and RvF model, we performed 10-fold cross validation to obtain a robust estimate of model performance. We report the mean area under the receiver operating curve (AUROC) across all of the k-folds as a measure of generalization performance effect size. To assess the statistical significance of the AUROC values, we used the p-value of the associated Mann-Whitney U statistic comparing the predicted scores between the true negative (trough/falling) and true positive (peak/rising) samples (**Figure S5,8**) [22]. We analyzed the magnitude and direction of the model coefficients to better understand the relationship between the EEG features and cycle phase. We applied the inverse PCA transform to the model coefficients to obtain individual feature importances. We tested for group-level trends in coefficient magnitude across the patient population using a two-tailed sign test against a zero median for each feature. #### General Statistical Tests For all paired (unpaired) non-parametric tests we used a sign test of differences (rank sum) to compare model effect sizes between patients. We used spearman correlation to analyze population trends between periodogram attributes and model performance. When comparing distributions within patients we applied parametric statistical tests for both paired and unpaired analyses. We tested for an uneven distribution of SEs at the population level both with and without stimulations across the four phase bins using a Kruskal-Wallis test, and for each individual we used a chi-square goodness of fit test against a uniform distribution to test for class imbalance. For all box plots, the middle line indicates the median, box edges are lower and upper quartiles, and whisker edges are the maximum and minimum. For all assessments, the significance threshold (α) was set at 0.05, and significance is indicated in figures with *= p < 0.05, **=p < 0.01, and \***| = p < 0.001. Because this paper is largely exploratory, we do not correct for multiple comparisons. ## Results ### Multidien cycles in interictal epileptiform activity In our retrospective cohort of 20 patients, we replicated analyses from previous studies [3], [5], [15] to compare the multidien dIEA cycle characteristics in our subjects with existing literature. We found significant multidien cycles in 89% of our patient population, with patient-specific cycle periods ranging from 7 to 49 days, consistent with prior studies [4], [5]. On average, subjects possessed 2 ± 0.65 significant peaks with a median significant multidien cycle length of 22.08 [13.13, 31.17] (**Figure S1**). Our analysis of the LE phase-locked cycles showed that the distribution of seizures across cycle phase was significantly non-uniform in 18 of 20 patients, with a median PLV of 0.24 [.17, .35], consistent with previously reported values (**Figure S2**)[4], [5]. We observed a significant negative correlation (R: -0.603, p = 0.005) between the peak period length of the LE phase-locked cycle and the LE phase-locking value, indicating stronger phase-locking in shorter multidien cycles (**Figure S3**). Because dIEA cycles represent fluctuations in stimulation event frequency, we quantified the disparity in the distribution of SE events both with *and* without stimulation across the four phase bins, given the higher number of stimulations at the cycle peak. Group-level analysis showed a significant difference in the distribution of available SE clips without stimulation across phase groups, with the largest difference between peak and trough recordings (KW test (72); Chi-sq = 33.8, p = 2.18e-07) (**Figure 1E**). When analyzing SEs containing stimulations, we found a smaller group-level class imbalance in the opposite direction (more peak than trough recordings) (KW test (67) Chi-sq = 8.53, p = 0.0363) (**Figure 1G**), (**Figure S4**). These results demonstrate distinct distributions between baseline SE recordings with and without stimulation. ### Spike rate is associated with dIEA counts We found that patients demonstrated a moderate correlation (Pearson *r* > 0.3) between spike rate and dIEA counts with a medium effect size (median *r*: 0.34, IQR: [.23, .48]). We next asked to what extent the spike rate extracted from the device recordings could differentiate between cycle phases. At the population level, we found that spike rate was greater at cycle peaks vs. troughs (sign test U(19)=15, p=0.0192), with 72% (15/19) of patients showing a higher spike rate at cycle peaks. Univariate effect sizes were not correlated with the positive predictive value of the spike detector (p > 0.5). We did not observe a significant difference at the population level in spike rate between the rising and falling dIEA cycle phases (**Figure 2E**). The spike rates we analyzed were calculated outside of SE recordings containing detected events, suggesting that the dIEA cycle is in part measuring the brain state that is associated with epileptiform spike rate in background recordings. ### Phase-dependence of interictal EEG features We asked whether features of background interictal activity were related to dIEA cycle phase (**Figure S5**). Our one-way ANOVA revealed that individual patients had at least one significant feature that was strongly associated with phase, though not the same feature for all patients. The proportion of day versus night recordings was similar across phase groups. The full results of the patient-specific ANOVA models are described in **Table S1 & Figure S7**. We next split our analysis into two paired comparisons, peak vs. trough and rising vs. falling, because we expected to see the largest differences in features between the peaks and the troughs of the cycles. For the PvT comparison, analysis of Cohen’s *d* effect sizes across patients revealed higher feature values at the peaks of dIEA cycles (**Figure 3A**). We observed significantly positive effect sizes for theta (sign test U(19) = 16, p = 0.004) and beta connectivity (U(19) = 17, p = 0.0007), as well as theta (U(19) = 16, p = 0.004) and alpha (U(19) = 16, p = 0.004) band power. Interestingly, we found a bimodal distribution of gamma band power effect sizes, indicating distinct patient subpopulations with higher gamma band power either at the peaks or troughs (**Figure 3B,S6**) Repeating the above analysis for RvF phases revealed that feature values were not preferentially higher in either phase, and the difference between rising and falling phases was not significant across the population in any feature (**Figure 3A, S7**). Our analysis of univariate features revealed no one electrophysiological biomarker that could strongly discriminate between the cycle peaks and troughs, but rather that subsets of features which discriminate between cycle phases are patient-specific. ![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/07/07/2023.07.05.23291521/F3.medium.gif) [Figure 3.](http://medrxiv.org/content/early/2023/07/07/2023.07.05.23291521/F3) Figure 3. Association of background features to dIEA cycle phase. A) Evolution of feature values over time for an example patient, trend lines are omitted for ≥5 days without recorded SE data (HUP096). Features were smoothed for visualization purposes and normalized between clinical visit (vertical lines). B) Heatmap demonstrating population trends in PvT effect sizes. The relationship between electrophysiology and dIEA is an individualized phenomenon, with different patients having different features that best associate with dIEA cycle phase. C) Population distributions of PvT and RvF multivariate model coefficients. PvT models reveal increased activity and connectivity during the peak compared to the troughs of multidien dIEA cycles. **D)** Mean out-of-sample PvT prediction ROC curve for HUP096 (purple) and ROC for each fold (light blue). **E)** Population distributions of model performance (AUROC) between PvT and RvF models without spike rate. Model performance is significantly higher in PvT models. ### ML Model with background features For each patient, we built binary linear SVM classifiers to predict cycle phase - PvT, and RvF - from the EEG features. Our models were able to predict PvT cycle phase above chance in all patients (AUC >0.5, 14/19 significant - **Figure S8, Table S1**), and RvF phase in 14/19 patients (8/19 significant - **Table S1**). Between the two classifiers, we saw higher performance in the PvT over the RvF model (sign test, p = 0.0044) (**Figure 3E**). Despite the fact that seizures preferentially occur in the rising phase of the cycle[5], our results indicate that EEG features are more strongly associated with the amplitude of dIEA cycles than individual phase bins. Even after accounting for correlated features (**Figure S9**), spike rate, and time of day, model coefficients showed higher band power and connectivity in the brain during the peaks of the multidien dIEA cycles (**Figure 3C**). While there was between-patient variation in the magnitude of the coefficients corresponding to each feature, we identified significant population trends of higher connectivity in the theta (sign test U(19) = 17, p = 0.0007) and beta (U(19) = 17, p = 0.0007) bands, and higher power in the theta (U(19) = 16, p = 0.0044) and alpha (U(19) = 17, p = 0.0007) bands at the peaks of the cycles. We also correlated model performance with both dIEA cycle length and LE phase-locking value, and found no significant trends. This suggests that the representation of dIEA cycles in interictal baseline electrophysiology is independent of cycle length as well as LE entrainment (**Figure S10**). ### Phase dependence of stimulation response features We hypothesized that multidien fluctuations would be temporarily superseded by the brain’s stimulation response, making multidien cycles more difficult to detect. Univariate analysis on post-stimulation periods was consistent with our univariate analysis of SE clips without stimulation, showing that brain connectivity and band power activity are increased at the peaks vs. troughs after stimulation (**Figure S11**). As with our background EEG analysis, we used an SVM model to explore the multivariate relationship between post-stimulation features and phase. The PvT detection model achieved above-chance performance in all patients (median AUC 0.58 [0.55, 0.64], **Figure 4B**) and the RvF model achieved above chance performance in 12 of 17 patients (median AUC 0.53 [0.50, 0.60], **Figure S12**), indicating that multidien phase remained detectable in the seconds following stimulation. However, compared to PvT detection in scheduled events without stimulation, post-stimulation detection performance was reduced, suggesting a relative suppression of the relationship between the EEG and dIEA cycles. Finally, the standardized feature importance for both PvT and RvF models did not exhibit a significant directional effect at the group level (**Figure 4A**). In general, both the direction and magnitude of feature importance is more heterogeneous across patients in the post-stimulation detection scenario, when compared with background EEG model coefficients. ![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/07/07/2023.07.05.23291521/F4.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2023/07/07/2023.07.05.23291521/F4) Figure 4: Post Stimulation Feature Importance and Model Comparison. **A)** Peak-Trough and Rising-Falling multivariate feature importance (Cohen’s D) in multivariate post-stimulation model. Stars indicate significance level of effect size distributions against a zero median (sign test). **B)** Comparison of PvT model performance (AUC) across SVM models trained with (1) background EEG features, (2) background EEG features + spikes, (3) spikes, (4) post-stimulation features. **C)** Comparison of EEG (no-stim) and post-stim model performance broken down by implant depth subgroups. Neo.: Neocortical implants, Mes.: Mesial temporal implants, Neo. + Mes.: one implant each in neocortical and mesial temporal structures. ### Comparison of ML model in different scenarios Finally, we compared PvT performance for models trained with the following feature sets: (1) features from SE clips without stimulation, (2) features from SE clips without stimulation AND spike rate, (3) spike rate, and (4) features from post-stimulation windows. Across the four feature sets, we found an increased ability to detect peak vs. troughs in the EEG (0.63 [0.56 - 0.67]) and EEG + spikes (0.62 [0.55 - 0.66]) models when compared to the spike only (0.55 [0.53-0.61], rank sum p = 0.0038) and post-stimulation models (0.58 [0.55-.64], rank sum p = 0.205) (**Figure 4B**). We further analyzed model performance by implant depth, the decrease in performance between background EEG to post-stim models appeared to be driven by patients with neocortical implants (median AUC 0.64 - 0.55) (**Figure 4C**). Similarly, the superior performance of PvT models when compared with RvF models persisted in the EEG feature sets but not in the post-stimulation or spike rate only feature set (sign test p= 0.03, 0.03, 1, 0.47). Although spike rate alone did not significantly indicate phase for all patients, we incorporated spike rate into the linear SVM of background features to determine whether a more complex interaction could better be associated with phase position. We found that the inclusion of spike rate did not substantially change the median AUC, though both the EEG and EEG with spikes models performed better than the spike rate only model, which had the worst performance across all models compared, indicating that the inclusion of background EEG features adds complementary information to phase detection (**Figure 4B**). ## Discussion Our study investigates RNS EEG recordings in 20 patients to build on prior work examining the relationship between EEG features and multidien cycles of cortical excitability. We leverage the unique recording and stimulation characteristics of the RNS dataset to complete the largest investigation to date on neuromodulation and dIEA cycles. Specifically, we show that: (1) beyond spike rate, features of interictal EEG signals contain information pertinent to detecting dIEA phase, (2) EEG features that are most associated with cycle phase are patient-specific, and (3) after neurostimulation, the EEG-dIEA relationship is suppressed. These findings may have important implications for tailoring pharmacological and neurostimulation therapy to patient chronotype, as well as for using multidien cycles as a predictive biomarker of seizure risk and understanding their underlying mechanisms. ### Spike rate does not tell the whole story: even after accounting for spikes, broader EEG features are associated with cycles Prior work investigating the multidien nature of epilepsy characterized the presence of cyclical seizure risk in acute event counts both in RNS dIEA[5], [12] as well as clinical spikes in different implantable devices[8]. We found evidence of multidien dIEA cycles in background EEG features beyond epileptiform spikes and acute signatures of epileptiform activity, and showed that multivariate EEG models outperformed spikes alone in predicting dIEA cycle phase. Our results suggest that there is complementary information between spikes and EEG features, and indicate that multidien cycles can be identified even below the spike detection threshold. It may be that some other latent mechanism, linked to both epileptiform activity and subthreshold EEG features, is responsible for cycle generation. EEG features could present an alternative or addition to using spikes in future therapeutic algorithms that incorporate patient chronotypes to titrate pharmacological or neuromodulatory therapy. ### Patient-specific features underlie the population level dIEA-EEG association The mechanism that causes variability in cycle period between patients is unknown and could be driving the observed differences in important feature sets between patients. It is also possible that differences in implant locations and electrode type could drive inter-patient variability in feature importance, given that measures of brain activity vary across brain structures [23]–[25]. Subgroup analysis revealed that the EEG-dIEA cycle coupling is strongest in neocortical patients. Prior work has shown that within patients dIEA cycles are stable between anatomical regions[5] but our finding suggests that there may be an anatomical basis for the presence of cycles in background recordings. Future work should rigorously investigate sub population effects on multidien cycles in background EEG features. ### The relationship between dIEA cycle and EEG features is attenuated after neurostimulation There is room for improving responsive stimulation specificity, and a functional question that arises from our observation of multidien cycles in epilepsy is: can stimulation perturb these cycles? Recently, Gregg et al. provided evidence that the power of the dominant dIEA cycle could be modulated by thalamic deep brain stimulation [7], and a broad body of literature has shown that stimulation acutely suppresses EEG signal features[17], [26], [27]. Our results are consistent with the potential for neurostimulation to disrupt cycles – we find a general suppression of the relationships between the EEG and dIEA with lower model effect sizes and limited population trends in model coefficients after stimulation, particularly with neocortical implants. One possible theory is that stimulation causes a stereotyped brain activity reset[27], reducing feature variance across cycle phases and model performance. Another theory is that the brain state triggering the stimulation is inherently decoupled from the dIEA cycles. ### Sampling bias in RNS recording types reflects cycle phase We found that cycle phase caused an imbalance in the availability of scheduled events both with and without stimulation, potentially biasing analyses towards recordings with different levels of cortical excitability. Given that values of EEG features fluctuate with cycle phase at the group level, sampling bias should be a practical consideration for studies separately analyzing RNS recordings with or without stimulation. ### Limitations Our study has several limitations. The relatively small size of our patient cohort, and the variability within it, precluded us from performing robust statistical sub-analyses or further exploring mechanistic hypotheses. Further studies with larger sample sizes are needed to confirm and expand upon our findings, and our group is currently working on a federated analysis of RNS data from multiple institutions to accomplish this goal [28]. We were also limited by the highly irregular sampling of EEG recordings by the RNS device, which prevented us from extracting cycles from EEG features using the same statistical tools and methods we applied to dIEA counts. This under sampling of the intracranial EEG is a fundamental challenge presented by this device. Future work investigating the observed phenomena should leverage continuous recording paradigms. Additionally, we select the cycle of interest based on LE PLV. While RNS LEs may not always capture true seizure activity, they have been shown to be correlated with electrographic events[29]. Despite these limitations, our study provides insight into the physiological correlates and underpinnings of the measured dIEA cycles in patients with RNS. #### Declaration of anonymous identifiers Sample/patient IDs disclosed in this study (HUPXXX, RNSXXX) are not identifiable and cannot be linked to any protected health information by any individuals outside of our research group. #### Declaration of generative AI and AI-assisted technologies in the writing process During the preparation of this work the authors used ChatGPT to edit for conciseness. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication. #### Declaration of competing interests E.C. is a paid consultant for Epiminder, an EEG device company, but declares no targeted compensation for this work. None of the other authors has any conflict of interest to disclose. ## Supporting information Supplemental Information [[supplements/291521_file03.docx]](pending:yes) ## Data Availability Data associated with the RNS device was obtained through a data use agreement between the University of Pennsylvania and NeuroPace, and cannot be publicly shared. However, code used to process RNS device data and perform our analysis can be found on GitHub at https://github.com/penn-cnt/RNS\_processing\_toolbox and https://github.com/wojemann/multidien-cycle-biomarkers [https://github.com/penn-cnt/RNS\_processing\_toolbox](https://github.com/penn-cnt/RNS_processing_toolbox) [https://github.com/wojemann/multidien-cycle-biomarkers](https://github.com/wojemann/multidien-cycle-biomarkers) ## Funding Statement This material is based on work was supported by the Mirowski Family Foundation, National Institute for Neurological Disorders and Stroke (NINDS) DP1NS122038, NINDS R61NS125568-01A1, NINDS K23NS121401, Pennsylvania Tobacco Fund, Burroughs Welcome Career Award for Medical Scientists, and National Science Foundation Graduate Research Fellowship under Grant No. DGE-1845298 ## Data Availability Statement Data associated with the RNS device was obtained through a data use agreement between the University of Pennsylvania and NeuroPace, and cannot be publicly shared. However, code used to process RNS device data and perform our analysis can be found on GitHub at [https://github.com/penn-cnt/RNS\_processing\_toolbox](https://github.com/penn-cnt/RNS_processing_toolbox) and [https://github.com/wojemann/multidien-cycle-biomarkers](https://github.com/wojemann/multidien-cycle-biomarkers) ## Acknowledgements We thank Magda Wernovsky for her help in consenting RNS patients. ## Abbreviations RNS : responsive neurostimulation dIEA : device-detected interictal epileptiform activity PLV : phase-locking value * Received July 5, 2023. * Revision received July 5, 2023. * Accepted July 7, 2023. * © 2023, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/) ## References 1. [1]. R. Bajorat, M. Wilde, T. Sellmann, T. Kirschstein, and R. Köhling, “Seizure frequency in pilocarpine-treated rats is independent of circadian rhythm,” Epilepsia, vol. 52, no. 9, pp. e118–e122, 2011, doi: 10.1111/j.1528-1167.2011.03200.x. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1528-1167.2011.03200.x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21801169&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000294973700005&link_type=ISI) 2. [2]. P. J. Karoly et al., “Cycles in epilepsy,” Nat. Rev. Neurol., vol. 17, no. 5, pp. 267–284, May 2021, doi: 10.1038/s41582-021-00464-1. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41582-021-00464-1&link_type=DOI) 3. [3]. M. G. Leguia et al., “Seizure Cycles in Focal Epilepsy,” JAMA Neurol., vol. 78, no. 4, p. 454, Apr. 2021, doi: 10.1001/jamaneurol.2020.5370. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/jamaneurol.2020.5370&link_type=DOI) 4. [4]. V. R. Rao, M. Leguia, T. K. Tcheng, and M. O. Baud, “Cues for seizure timing,” Epilepsia, vol. 62, no. S1, Feb. 2021, doi: 10.1111/epi.16611. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/epi.16611&link_type=DOI) 5. [5]. M. O. Baud et al., “Multi-day rhythms modulate seizure risk in epilepsy,” Nat. Commun., vol. 9, no. 1, p. 88, Dec. 2018, doi: 10.1038/s41467-017-02577-y. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-017-02577-y&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29311566&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) 6. [6]. P. J. Karoly et al., “Circadian and circaseptan rhythms in human epilepsy: a retrospective cohort study,” Lancet Neurol., vol. 17, no. 11, pp. 977–985, Nov. 2018, doi: 10.1016/S1474-4422(18)30274-6. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1474-4422(18)30274-6&link_type=DOI) 7. [7]. N. M. Gregg et al., “Thalamic deep brain stimulation modulates cycles of seizure risk in epilepsy.,” Sci. Rep., vol. 11, no. 1, Dec. 2021, doi: 10.1038/s41598-021-03555-7. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41598-021-03555-7&link_type=DOI) 8. [8]. M. I. Maturana et al., “Critical slowing down as a biomarker for seizure susceptibility,” Nat. Commun., vol. 11, no. 1, p. 2172, May 2020, doi: 10.1038/s41467-020-15908-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-020-15908-3&link_type=DOI) 9. [9]. N. M. Gregg et al., “Circadian and multiday seizure periodicities, and seizure clusters in canine epilepsy,” vol. 2, no. 1, Feb. 2020, doi: 10.1093/braincomms/fcaa008. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/braincomms/fcaa008&link_type=DOI) 10. [10]. M. O. Baud, A. Ghestem, J.-J. Benoliel, C. Becker, and C. Bernard, “Endogenous multidien rhythm of epilepsy in rats,” Exp. Neurol., vol. 315, pp. 82–87, May 2019, doi: 10.1016/j.expneurol.2019.02.006. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.expneurol.2019.02.006&link_type=DOI) 11. [11]. P. J. Karoly et al., “Multiday cycles of heart rate modulate seizure likelihood at daily, weekly and monthly timescales: An observational cohort study,” Neurology, preprint, Nov. 2020. doi: 10.1101/2020.11.24.20237990. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoibWVkcnhpdiI7czo1OiJyZXNpZCI7czoyMToiMjAyMC4xMS4yNC4yMDIzNzk5MHYyIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDcvMDcvMjAyMy4wNy4wNS4yMzI5MTUyMS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 12. [12]. T. Proix et al., “Forecasting seizure risk in adults with focal epilepsy: a development and validation study,” Lancet Neurol., vol. 20, no. 2, pp. 127–135, Feb. 2021, doi: 10.1016/S1474-4422(20)30396-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1474-4422(20)30396-3&link_type=DOI) 13. [13]. P. J. Karoly et al., “Circadian and circaseptan rhythms in human epilepsy: a retrospective cohort study,” Lancet Neurol., vol. 17, no. 11, pp. 977–985, Nov. 2018, doi: 10.1016/S1474-4422(18)30274-6. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1474-4422(18)30274-6&link_type=DOI) 14. [14]. M. J. Morrell, “Responsive cortical stimulation for the treatment of medically intractable partial epilepsy,” Neurology, vol. 77, no. 13, pp. 1295–1304, 2011, doi: 10.1212/WNL.0b013e3182302056. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1212/WNL.0b013e3182302056&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21917777&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) 15. [15]. D. N. Anderson et al., “Closed-loop neurostimulation for epilepsy leads to improved outcomes when stimulation episodes are delivered during periods with less epileptiform activity,” Neurology, preprint, Nov. 2022. doi: 10.1101/2022.11.28.22282784. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoibWVkcnhpdiI7czo1OiJyZXNpZCI7czoyMToiMjAyMi4xMS4yOC4yMjI4Mjc4NHYxIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDcvMDcvMjAyMy4wNy4wNS4yMzI5MTUyMS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 16. [16]. M. W. Brown et al., “Comparison of novel computer detectors and human performance for spike detection in intracranial EEG,” Clin. Neurophysiol., vol. 118, no. 8, pp. 1744–1752, Aug. 2007, doi: 10.1016/j.clinph.2007.04.017. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.clinph.2007.04.017&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17544322&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000248668600012&link_type=ISI) 17. [17]. S. N. Rønborg et al., “Acute effects of brain-responsive neurostimulation in drug-resistant partial onset epilepsy,” Clin. Neurophysiol., vol. 132, no. 6, pp. 1209–1220, Jun. 2021, doi: 10.1016/j.clinph.2021.03.013. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/J.CLINPH.2021.03.013&link_type=DOI) 18. [18]. P. J. Karoly et al., “Multiday cycles of heart rate are associated with seizure likelihood: An observational cohort study,” eBioMedicine, vol. 72, p. 103619, Oct. 2021, doi: 10.1016/j.ebiom.2021.103619. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ebiom.2021.103619&link_type=DOI) 19. [19]. C. Torrence and G. P. Compo, “A practical guide to wavelet analysis,” Bull. Am. Meteorol. Soc., vol. 79, no. 1, pp. 61–78, Jan. 1998, doi: 10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2&link_type=DOI) 20. [20]. P. Berens, “CircStat: A MATLAB Toolbox for Circular Statistics,” J. Stat. Softw., vol. 31, pp. 1–21, Sep. 2009, doi: 10.18637/jss.v031.i10. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.18637/jss.v031.i10&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25402854&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) 21. [21]. M. Lin, H. C. Lucas, and G. Shmueli, “**Research Commentary** —Too Big to Fail: Large Samples and the *p* -Value Problem,” Inf. Syst. Res., vol. 24, no. 4, pp. 906–917, Dec. 2013, doi: 10.1287/isre.2013.0480. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1287/isre.2013.0480&link_type=DOI) 22. [22]. S. J. Mason and N. E. Graham, “Areas beneath the relative operating characteristics (ROC) and relative operating levels (ROL) curves: Statistical significance and interpretation,” Q. J. R. Meteorol. Soc., vol. 128, no. 584, pp. 2145–2166, 2002, doi: 10.1256/003590002320603584. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1256/003590002320603584&link_type=DOI) 23. [23]. J. M. Bernabei et al., “Normative intracranial EEG maps epileptogenic tissues in focal epilepsy,” Brain, vol. 145, no. 6, pp. 1949–1961, Jun. 2022, doi: 10.1093/brain/awab480. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/brain/awab480&link_type=DOI) 24. [24]. A. N. Khambhati, A. Shafi, V. R. Rao, and E. F. Chang, “Long-term brain network reorganization predicts responsive neurostimulation outcomes for focal epilepsy,” Sci. Transl. Med., vol. 13, no. 608, p. eabf6588, Aug. 2021, doi: 10.1126/scitranslmed.abf6588. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTE6InNjaXRyYW5zbWVkIjtzOjU6InJlc2lkIjtzOjE1OiIxMy82MDgvZWFiZjY1ODgiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMy8wNy8wNy8yMDIzLjA3LjA1LjIzMjkxNTIxLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 25. [25]. P. N. Taylor et al., “Normative brain mapping of interictal intracranial EEG to localize epileptogenic tissue,” Brain, vol. 145, no. 3, pp. 939–949, Mar. 2022, doi: 10.1093/brain/awab380. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/brain/awab380&link_type=DOI) 26. [26]. A. N. Khambhati et al., “Functional control of electrophysiological network architecture using direct neurostimulation in humans,” Netw. Neurosci., vol. 3, no. 3, pp. 848–877, Jan. 2019, doi: 10.1162/netn\_a\_00089. [CrossRef](http://medrxiv.org/lookup/external-ref?access\_num=10.1162/netn_a_00089&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31410383&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) 27. [27]. V. S. Sohal and F. T. Sun, “Responsive Neurostimulation Suppresses Synchronized Cortical Rhythms in Patients with Epilepsy,” Neurosurg. Clin. N. Am., vol. 22, no. 4, pp. 481–488, Oct. 2011, doi: 10.1016/j.nec.2011.07.007. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.nec.2011.07.007&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21939847&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F07%2F07%2F2023.07.05.23291521.atom) 28. [28]. B. H. Scheid et al., “Intracranial electroencephalographic biomarker predicts effective responsive neurostimulation for epilepsy prior to treatment,” *Epilepsia*, p. epi.17163, Jan. 2022, doi: 10.1111/epi.17163. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/epi.17163&link_type=DOI) 29. [29]. M. Quigg, T. L. Skarpaas, D. C. Spencer, N. B. Fountain, B. Jarosiewicz, and M. J. Morrell, “Electrocorticographic events from long-term ambulatory brain recordings can potentially supplement seizure diaries,” Epilepsy Res., vol. 161, no. October 2019, p. 106302, Mar. 2020, doi: 10.1016/j.eplepsyres.2020.106302. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.eplepsyres.2020.106302&link_type=DOI) 30. [30]. M. Farge, “Wavelet transforms and their applications to turbulence,” Annu. Rev. Fluid Mech., vol. 24, pp. 359–457, 1992. 31. [31]. W. F. Massy, “Principal Components Regression in Exploratory Statistical Research,” J. Am. Stat. Assoc., vol. 60, no. 309, pp. 234–256, Mar. 1965, doi: 10.1080/01621459.1965.10480787. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1080/01621459.1965.10480787&link_type=DOI)