Seizure risk in epilepsy changes with multi-day biological rhythms

We recruited 37 subjects (22 males) who had been implanted with the RNS System for purely clinical indications for at least 3 months and up to 9.9 years (median: 2.25 years) across two neurology centers (University of California, San Francisco, N = 11, and California Pacific Medical Center, N = 26). The two Institutional Review Boards approved the study and written informed consent was obtained from all subjects. Subjects had a variety of focal epilepsies (Supplementary Table 1). Indications for treatment with the RNS System as opposed to resective surgery included bilateral seizure localization (temporal and frontal), seizures arising from eloquent cortex (motor and visual), and previous contralateral resection.

Detection of epileptiform activity by the RNS System relies on user-configurable tools (line length, area under the curve, and bandpass filtering) for which thresholds are optimized to detect seizure onset patterns⁴². Examples of subject-specific epileptiform activity detected by the RNS System have been reported previously³ and are shown in Supplementary Fig. 1. Hourly detection counts are stored by the RNS System for the last 28 days and are continuous as long as subjects download device data at least this often. For each subject, we discarded data recorded from the day of implant until reliable detection of seizures and IEA was achieved during the first few outpatient visits (median number of visits: 2, range: 1–11) a few months later (median: 29 days, range: 2–90 days; Supplementary Fig. 4). Sixteen subjects who did not download device data regularly had resultant gaps in detection counts. The data containing gaps longer than 6 days was considered discontinuous and analyzed in separate segments. The data of <90 days surrounded by gaps was discarded. Gaps up to 6 days were interpolated (see below). One subject (S4) underwent resective surgery and the device continued to record. This data were discarded. For visualization purposes daily counts were obtained by averaging 24 h of data on the same calendar date, but actual analyses were performed on the original hourly count data.

Due to memory constraints, the RNS System can store only a limited number of raw EEGs at a given time. Seizure detection with the RNS System relies instead on a surrogate marker, long-epileptiform activity (LEA; also termed ‘Long Episode’ by others^3,42), which occurs when the EEG signal meets detection criteria for a clinician-specified length of time (typically, the minimum duration of each subject’s electrographic seizure; LEA durations for this study ranged from 15 to 40 s, with an average of 26.8 s). Like detection counts, LEA timestamp information is stored by the device and available for analysis. Although generally a reliable proxy for seizures⁴³, LEA can also represent epochs of abundant IEA that do not meet criteria for electrographic seizures⁴⁴ (Supplementary Fig. 1b). To avoid contaminating seizure analysis with such false positives, we calculated for each subject the positive predictive value (PPV) of LEA for electrographic seizures (Supplementary Fig. 1c). For each subject, a Board-certified epileptologist (V.R.R.) assessed whether individual LEA corresponded to true electrographic seizures by visually reviewing corresponding EEGs. When possible (<300 LEA, 10 subjects), all individual EEGs were reviewed and each was labeled as an electrographic seizure (true positive) or as LEA other than an electrographic seizure (false positive). For subjects with too many LEA EEGs for comprehensive review (N > 300, 27 subjects), 100 EEGs were randomly selected from epochs with stable detection settings and a PPV was calculated for each epoch. Subjects (N = 3) who had <20 LEA were excluded from statistical analysis, because this was too low a number for a histogram-based statistical technique (see below). In total, we included 14 subjects for whom LEA was a reliable surrogate for seizures (PPV >90%; mean 98%, range 92–100%; Supplementary Table 1). Thus, we estimate that <2% of LEA used here may actually be sustained trains of IEA and that the rest represent true electrographic seizures.

Hourly IEA counts were normalized (z-score) separately by block of recording between clinic visits so as to ensure stability of detection settings and anticonvulsant medications for each block. Continuous time-series were obtained by concatenating these individual blocks. Power and phase of the hourly IEA counts time-series were obtained using a Morlet wavelet transform for 89 period bins (scales) with increasing spacing: 1.2 h between 2.4 and 31.2 h, 2.4 h between 33.6 and 48 h, 4.8 h between 2.2 and 4 d, 12 h between 4.5 and 10 d, and 24 h between 11 and 45 d. Gaps in recordings, if relatively short, were interpolated by a method similar to kriging in geostatistical and climatology research⁴⁵ (Supplementary Fig. 3). For each gap, the variance was calculated for two windows, before and after the gap (each with same length as the gap itself). This variance was used to cast the normalization curves for pseudo-random data generated around a central tendency line directly connecting the means of the two peri-gap windows above. This process was performed for each gap, and only for gaps shorter than 20% of the period being analyzed (e.g., up to 2 days of interpolation allowed for a period of 10 days). Gaps were processed in increasing duration order, so that higher frequencies were not affected by interpolations performed at lower frequencies. In addition, a cone of influence—the region of the wavelet spectrum, shaped according to period length, in which edge effects impede accurate periodic estimation—was discarded at the extremities and around gaps too wide to be interpolated. The data could then be represented as a spectrogram of power or phase over time (time-frequency analysis) with excluded data at the extremities and around gaps (Supplementary Fig. 2). Power index in periodograms was estimated for each individual scale (period bin) as the square root of the average over time of the absolute value of complex wavelet coefficients. Principal component analysis of all individual periodograms (Supplementary Fig. 3) was used to extract recurrent patterns. The six first components were selected as they explained 98.7% of the total variance in the frequency domain. Average periodograms were calculated for three separate clusters obtained by K-means (cosine distance, rank of 3) on the principal components coefficients (Fig. 2). Peaks in periodicity were defined as a positive-to-negative zero-crossing of the derivative of the periodogram.

Phase at every time-point was calculated in the frequency domain for a single band of wavelet coefficients corresponding to the peak period ±33.3% (e.g., 24 ± 8 h, or 15 ± 5 d), so as to accommodate variation in periodicity. Thus, at each time-point, the most powerful frequency within a range would most influence the phase. The average and variance in period length were calculated using the distance between two successive phase values of zero, excluding gaps. Circadian epileptiform peak activity histograms were obtained by counting occurrences within 1-h bins (24 bins, Fig. 3). Seizure phase histograms were obtained by counting occurrences within 20-degree phase bins (18 bins from –π to +π) for circadian and multidien rhythms (Fig. 4). Similarly, average amplitude was obtained by calculating the mean and standard deviation of normalized data for all time-points with phase contained within these bins. If more than one multidien peak was present in the periodogram, we used the shortest (first) multidien periodicity to evaluate instantaneous phase. For visualization purposes, signal was reconstructed for the peak periods ±33.3% using an inverse wavelet transform (Fig. 1 and Supplementary Fig. 2). The phase-space representation allowed for the study of seizure risk as it relates to circadian vs. multidien phase. Scatterplots help visualize these relations abstracted from the fact that multidien periods vary across subjects.

The sample size was determined by the availability of data. Previous studies have demonstrated quantitative analyses of these data sets with a relatively small number of subjects⁵ (though our sample size is considerably larger in both number of subjects and recording durations). Circular statistics can be applied within-subject, and we investigated replicability across all subjects in our data set. The only pre-established exclusion criterion was length of continuous recording of <3 months, and we excluded seven recently implanted subjects for this reason. Given the observational nature of the study, there was no replication of measurements over time per se, though these long recordings and autocorrelation analyses (Supplementary Fig. 5) serve as technical measures of within-subject replicability. Values were expressed as mean ± standard deviation (SD) and plotted as dots and error-bars, unless specified otherwise. For categorical distributions (gender and seizure localization), we performed a χ²-test. For continuous variables, we performed a Wilcoxon test. Autocorrelations to the average periodogram were assessed at each time-point in the frequency domain using a Pearson coefficient and averaged over the total length of recording for each individual. For eight subjects with bilateral implants and >3 months of bilateral recordings, Pearson coefficients were calculated on the periodograms derived from each site. Other correlations were assessed using a linear model with intercept and the F-statistic vs. a constant model. Seizure rate was log-transformed for regression models, as it had a logarithmic distribution. Phase analyses were done using the 2012 Matlab circular statistics toolbox by Dr. Philipp Berens⁴⁶ including functions for circular mean, circular standard deviation, and resultant vector. We used the circular variant of k-means for clustering of angles. We used the Omnibus (or Hodges-Ajne) test to calculate statistical significance for non-uniform angular distribution (against the null hypothesis of a uniform distribution), as opposed to the more classical Rayleigh test, because some angular distributions were bimodal, especially in circadian rhythms. We used the Kuiper test (analog to the Kolmogorov–Smirnov test for circular data) to calculate the statistical difference between two angle populations, as it makes no assumption on the underlying distributions. These descriptive statistics, clustering methods, and tests in the circular domain have equivalents in the linear domain but they take into account the fact that angles flip from π to –π, creating numerical discontinuity. The mean resultant vector is a metric of phase consistency of a population of angles. It is constructed by vectorial averaging of phase vectors each of unitary length equal to one. Thus, the angle of the mean resultant vector represents the mean angle of the population and its length (Phase-Locking Value, PLV) represents dispersion (values close to zero) or concentration (values close to one) of constituting angles. Given that the PLV is a continuous variable linearly varying from zero to one, we used a Wilcoxon test to compare population medians. To study the stability of this metric over time, we computed repeated PLVs every 3 months including data over a year (±6 months at each time-point, Fig. 4f). To study the map of seizure risk for each individual subject, we binned seizure counts in the phase-space into 324 (18 × 18) 20-degrees circadian and multidien phase combination elements. To reduce noise, we smoothed this result with a 2D Gaussian kernel spanning two standard-deviations over ±40 degrees. Each single bin, containing the number of seizures that occurred at a given phase combination, successively represented true positives (TP), whereas seizure counts in the 323 other bins represented false negatives (FN). False positive (FP) were the number of occurrence of a given phase combination (without seizures), and True Negative (TN) were all the other phase combinations (without seizures). The risk ratio was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\rm TP}/({\rm TP + FP})}}{{{\rm FN}/({\rm TN + FN})}}$$\end{document}TP∕(TP+FP)FN∕(TN+FN). This was done for each bin resulting in a map of ‘risk factor-seizure’ association (Fig. 5). Thus, the risk ratio takes into account the number of seizures observed at a given phase combination, the seizure rate of a given subject, and the probability of finding a given combination of phases. These maps were averaged across subjects after alignment to the preferred circadian and multidien phases (Fig. 6). In addition to this visualization of risk-ratio maps, we estimated the global effect size with confidence intervals for each individual, by grouping bins of the map that were in-phase (the ‘risk factor’) or anti-phase with the preferred angle (±90°) and calculating a risk-ratio between these two categories (Fig. 5). This was done for circadian and multidien rhythms independently and for the combination of the two phases (using a union Boolean operator). Finally, we displayed risk ratios as a Forest plot and computed an effect size summary metric across subjects using two methods: (i) a simple unweighted average of the risk ratio and standard error, and (ii) a random effect weighted average of the risk ratio and standard error (variance in risk-ratio was assumed to be non-random, Fig. 5)⁴⁷.

The data and code utilized in this study are available from the corresponding author upon reasonable request.

Supplementary Information Peer review File