Fluctuations in EEG band power at subject-specific timescales over minutes to days are associated with changes in seizure dynamics
Mariella Panagiotopoulou, Christoforos Papasavvas, Gabrielle M Schroeder, Peter Taylor, Yujiang Wang
FFluctuations in EEG band power over minutes to days explain how seizures change over time
Mariella Panagiotopoulou , Christoforos Papasavvas , Gabrielle M Schroeder ,Peter Taylor , , , Yujiang Wang ∗ , , December 15, 2020
Abstract
Epilepsy is recognised as a dynamic disease, where seizures and their features change overtime. Specifically, we recently demonstrated that seizures themselves change in terms of theirevolution. However, the underpinning modulators of seizure evolution remain unclear.In this work, we analyse continuous (interictal) intracranial Electroencephalographic (iEEG)recordings, and elucidate fluctuations in iEEG band power over different timescales (rangingfrom minutes to days).We find that all patients show an approximately circadian fluctuation in their EEG bandpower, but also many other fluctuations on patient-specific timescales. Importantly, we findthat a combination of fluctuations on different timescales can explain the changes in seizureevolution in most patients above chance-level.We interpret these results as evidence that seizure modulating factors exist, and they varyover time (patient-specifically). These time-varying modulating factors can be captured influctuations of EEG band power, and future work should link them to the exact biologicaltime-varying processes. a r X i v : . [ q - b i o . N C ] D ec Introduction
Epilepsy is a common neurological condition, characterised by recurrent, unprovoked seizures [1]. Itaffects approximately 1% of the world’s population and a third of the patients experience refractoryepilepsy, where seizures are not adequately controlled despite medication.Importantly, epilepsy is not a static disorder; seizures and pathological brain activity have beenshown to fluctuate over hours to years. In terms of the seizures, changes are observed in terms oftheir evolution, spread and dynamics from one seizure to the next. Clinical symptoms and severityof the seizure also correlate with the observed changes in seizure dynamics. While a subset ofseizures might share common features in the same patient, our recent work [2] has shown thatfunctional network dynamics of epileptic seizures change over time in the same patient. Notably,these changes followed daily (circadian) fluctuations in approximately a third of the patients, andlonger-term fluctuations, or a mixture between the two was seen in almost all remaining patients[2]. In support of our observations, a recent study quantifying univariate properties of seizureonset and offset also noted that different types of dynamics can be seen across different seizuresin the same patient [3, 4]. Thus, epileptic seizures are not a deterministic sequence of abnormalbrain activity patterns, but are clearly modulated by processes, which shape the neural activityduring a seizure and affect the severity of the seizure.One hypothesis explaining these observations is that the diversity of seizures may arise fromfluctuations in the underlying interictal brain states. It is well-known that spectral properties ofthe electroencephalogram (EEG) change in each channel from moment to moment [5]. Global andlocal characteristics of the interictal functional network have been demonstrated to fluctuate overthe timescale of hours to days, with a strong effect of the circadian rhythm [6, 7, 8]. Interictalfluctuations related to epilepsy have also been reported. For example, high frequency oscillation(HFO) rates vary in location and power within each patient over time [9]. Interictal spikes alsochange in their location and rate over hours to days [10, 11]. Interestingly, periodicity in interictalspikes over long time scales (days, weeks, months) has also been shown to be associated with seizureoccurrence [10, 12]. Therefore, it remains an open question if variability in seizure dynamics canbe explained to some degree by temporal fluctuations of interictal brain activity on different time-scales.We seek to address this question by exploring what temporal fluctuations of interictal brainactivity best explain the diversity in seizure dynamics for each individual patient. Most previousresearch has focused on analysing spikes, bursts, or HFOs and their association with seizure oc-currence. However, an analysis of the factors driving changes in dynamics from one seizure to thenext is missing. Additionally, we seek to explore the full range of all possible brain activity pat-terns, which will include epilepsy-related activity patterns (e.g. spikes), but are patient-specific.Therefore, we will analyse temporal fluctuations in wide-ranging spectral properties of long-termEEG, using with patient-specific pattern detection. We will analyse which frequency bands, orgroups of channels contribute to fluctuations of particular timescales. Finally, we will explore theassociation of fluctuations on different timescales with seizure variability in each patient.
For the purpose of this study, we used open source data from patients with drug-resistant epilepsy(available at http://ieeg-swez.ethz.ch). The data consist of a total of 2656 hours of long-term in-2 atient Duration Duration Number of Mean seizurein hours in days Seizures duration (minutes)ID01
293 12 2 10.030
ID02
235 10 2 1.468
ID03
158 7 4 1.078
ID04
41 2 14 0.699
ID05
109 5 4 0.278
ID06
146 6 8 0.765
ID07
69 3 4 1.159
ID08
144 6 70 0.366
ID09
41 2 27 0.706
ID10
42 2 17 1.181
ID11
212 9 2 1.526
ID12
191 8 9 2.441
ID13
104 4 7 1.717
ID14
161 7 60 0.430
ID15
196 8 2 1.576
ID16
177 7 5 3.174
ID17
130 5 2 1.632
ID18
205 9 5 3.319
Total
Average
Table 1:
For each subject (
Patient ) the following descriptive information is provided:
Durationin hours: duration of iEEG recordings in hours.
Duration in days: duration of iEEG recordingsin days.
Number of Seizures: number of patient’s seizures annotated.
Mean seizure duration: mean duration across all seizures in minutes.tracranial electroencephalography (iEEG) from 18 subjects. Continuous recordings in each subjectcover 24 to 128 channels and vary between 2 to 12 days. Sampling frequency was either 512 or1024 Hz depending on the subject. Electrodes (strip, grid, and depth) were implanted intracra-nially by clinicians. The collection of the data was conducted in the Sleep-Wake-Epilepsy-Center(SWEC) at the University Hospital of Bern, Department of Neurology, as part of their presurgicalevaluation programme, independently of this study [13].IEEG signals were provided in already pre-processed form. Briefly, signals have been median-referenced and band-pass filtered from 0.5-120 Hz using a th order Butterworth filter (forwardand backward). Seizure onset and termination times were determined by a board-certified epilep-tologist. Channels with artifacts were also identified and excluded by the same epileptologist. Allthe steps so far were conducted independently of this study resulting to publicly available data.All patients formally consented to their iEEG data being used for research purposes[13]. We performed additional preprocessing steps to extract iEEG band power from the five mainfrequency bands (Fig. 1a). For each recording channel, the signal was divided into 3 s epochs(Fig. 1b). For each epoch, the band power was computed for the following frequency bands:3 a) t i m e (c) (d) NMF Decomposition Time (days) C h a nn e l a nd f r e qu e n cy b a nd (b) iEEG Signals Band power over time ... ch1ch32 d e l t a t h e t aa l ph a b e t a g a mm a l og band po w e r , s t anda r d i s ed , s i g m o i d t r an s f o r m ed ComponentsW - component weights de l t a t he t aa l phabe t aga mm a ch1...ch32 H - component expression coefficientsTime (days) C o m pon e n t s Figure 1: Workflow of data pre-processing; calculation of Band Power and NMFimplementation . (a) 30 s segments were extracted along the recording from the multi-channeliEEG. The presentation of iEEG signals here is used as an example; it does not reflect 30 s recordingsegments. (b) The log of the standardised band power. (c)&(d) NMF is performed to the bandpower matrix resulting into the decomposition W × H . (c) W contains the basis vectors; (d) Thecoefficient matrix H gives the pattern of the contribution of each basis vector to the observed valuesof frequency band. δ : 1 − Hz, θ : 4 − Hz, α : 8 − Hz, β : 13 − Hz and γ : 30 − Hz. The band powers werecomputed using Welch’s method with 1 s sliding window and an overlap of 0.5 s between consecutivewindows. To ensure computational tractability, the band power was subsequently averages every30 s (i.e. every 10 epochs) giving rise to a time series of band power that was sampling an averageevery 30 s. The band power values were aggregated into band-specific matrices with dimensions4channels × S ( x ) = [1 + exp( − x )] − ) the standardised data to ensure positive entries between 0and 1. For each patient, we then concatenated the matrices from all frequency bands yielding asingle (5 × × n × T ) matrix (Fig. 1b). We will alsorefer to this matrix as the data matrix X throughout the paper. As the data matrix X is high-dimensional with redundant information, we applied a dimensionalityreduction step to the data. Dimensionality reduction was performed on the matrix X defined inSection 2.2 using NMF [14]. NMF approximates a non-negative input matrix X ∈ R n × T + as theproduct of two non-negative matrices, W ∈ R n × k + and H ∈ R k × T + , such that X ≈ W × H ≡ X (cid:48) ,given an integer k . Specifically, we applied the non-negative singular value decomposition (SVD)with low-rank correction (NNSVD-LRC) [15], which is an NMF algorithm based on SVD that usesin every iteration of the optimisation process the discarded SVD factors for the initialisation step.In this way, we can decompose each subject’s band power data matrix X into W and H matrices(Fig. 1c-d). Every column of matrix W represents a single NMF component and is a basis vector orfeature vector with n elements/features. Each row of H represents how a single NMF componentevolves over time across all T time epochs.To determine the optimal number of representative NMF components, k , we scanned ≥ k ≥ in each subject. For each value of k we performed NNSVD-LRC, obtained a W and H and calculated the relative reconstruction errors (cid:80) n,T | X − X (cid:48) | / ( n × T ) . We also calculated c = max { max ( | Corr ( W ) | ) , max ( | Corr ( H ) | ) } for each k , where max ( | Corr ( W ) | ) represents themax absolute correlation among all column pairs of W , and max ( | Corr ( H ) | represents the maxabsolute correlation among all row pairs of H . In simple words, we calculated the strongestcorrelation or anticorrelation between NMF components in terms of their feature weights W andtheir expression coefficient time series H . We then selected k = 15 , or the k yielding the smallestcorrelation between the NMF components that has a relative reconstruction error smaller than , which ever is smaller. This means that different subjects had a distinct number of NMFcomponents, k . After applying NMF with the optimal choice of k , we obtained two matrices for each subject, W and H . The matrix W consists of the basis vectors, while H is a multivariate time-series withdimensions equal to k × T ( = the number of NMF components × the total number of epochs), asalready described in Section 2.3.To investigate fluctuations in the band power data on different timescales, we analysed thematrix H using Empirical Mode Decomposition (EMD) [16, 17]. It is well known that EEGsignals are non-stationary processes characterised by time-varying features [18, 19]. EMD is adata-adaptive method able to detect non-stationary and non-rhythmic fluctuations on differenttimescales. EMD decomposes an input signal, say Y ( t ) , into M finite narrow-band fluctuations,known as intrinsic mode functions (IMF), based on the local extrema of a time-series: Y ( t ) = (cid:80) Mi =1 IM F i ( t ) + r ( t ) , where r ( t ) is the residue signal. This process is called sifting process [16].For the purpose of this analysis, we used an extension of this method to the multi-dimensionalspace called Multivariate Empirical Mode Decomposition (MEMD) [16, 17, 20]. Therefore, the5atrix H can be represented by the sum of M multi-dimensional IMF signals (dimension foreach IMF would be equal to k ; i.e. the number of rows of the matrix H /NMF components). Toclarify, we can think of all IMFs in a specific dimension as the decomposition of the correspondingNMF component. If we denote IM F i,j the dimension j of the IM F i signal, then the j -th NMFcomponent Y j ( t ) can be written as Y j ( t ) = (cid:80) Mi =1 IM F i,j ( t ) + r j ( t ) . This equation applies to everyNMF component j = 1 , . . . , k .We then apply a Hilbert-transform on each dimension of the IMF (following classical analysismethods for EMD) to obtain the instantaneous frequency, phase, and amplitude.For any (real-valued) univariate signal, say u ( t ) , we can define its Hilbert transform as: H ( u )( t ) = 1 π P (cid:90) + ∞−∞ u ( τ ) t − τ dτ, (1)where P represent the Cauchy principal value for any function u ( t ) ∈ L P class [16].The analytical signal v ( t ) obtained from the Hilbert transform can be expressed as: v ( t ) = u ( t ) + iH ( u )( t ) = a ( t ) e iθ ( t ) , (2)where a ( t ) = (cid:112) u ( t ) + H ( u )( t ) (3)and θ ( t ) = tan − (cid:18) H ( u )( t ) u ( t ) (cid:19) (4)where a ( t ) is the instantaneous amplitude and θ ( t ) is the instantaneous phase, respectively.The instantaneous frequency, f ( t ) can then be calculated as follows: f ( t ) = dθ ( t ) dt . (5)The application of EMD along with Hilbert transform leads to the so-called Hilbert-Huangtransform. Through the Hilbert spectral analysis, each IMF’s instantaneous frequency can berepresented as functions of time. The result is a frequency-time distribution of signal amplitude(or energy using the squared values of amplitude, a ( t ) ), designated as Hilbert amplitude spectrumor Hilbert spectrum (or Hilbert energy spectrum if energy is used instead of amplitude), H ( f, t ) .For each univariate IMF, we can obtain the Hilbert energy spectrum as a function of instanta-neous frequency and time mathematically using the following formulas: H ( f, t ) = (cid:26) a ( t ) , f = f ( t )0 , otherwise. (6)The Hilbert-Huang marginal spectrum h ( f ) of the original signal Y ( t ) can then be definedas the total energy distributed across the frequency space within a defined time period [0 , T ] .Mathematically, this can be expressed as shown below: h ( f ) = (cid:90) T H ( f, t ) dt. (7)By using Equations 6 & 7 we can obtain the Hilbert-Huang marginal spectrum for a univariateIMF. However, the application of the multivariate EMD results in multivariate IMF signals. Inorder to compute the Hilbert-Huang marginal spectrum of each IMF across all dimensions, wesimply average H i ( f, t ) across i = 1 , . . . , k dimensions.6 .5 Relative power in each IMF To understand how much each of the iEEG frequency bands contributed to a certain IMF, we firstdetermined how much each dimension of the IMF contributed to the overall power of the IMF. Tothis end, we first obtained the mean power E ij in each dimension j of every i -th IMF signal: E ij = (cid:80) Tt =0 a ij ( t ) T , (8)where T is, as before, the number of time epochs, and a ij ( t ) is the instantaneous amplitude for the j -th dimension of i -th IMF signal at time point t . One of the main properties of MEMD is thatmultivariate signals are decomposed in multivariate IMF signal of the same dimensions, where alldimensions within an IMF share the same time-scale of fluctuation [21]. Hence, focusing on themean power over time of each dimension within an IMF is a good indication of the power on aparticular time scale. The relative contribution of each j -th dimension to the i -th IMF (or relativepower) is then simply R ij = E ij (cid:80) kj =1 E ij , (9)with k indicating the number of dimensions. In order to check if all channels contribute homogeneously to a NMF component in a particularfrequency band, we used a measure of sparsity: the
Gini Index or Gini coefficient [22]. Given avector x = ( x , x , . . . , x N ) sorted in ascending order such that x < x < · · · < x N , the Gini Indexcan be derived using the following formula: G ( x ) = 1 − N (cid:88) i =1 x i (cid:107) x (cid:107) (cid:18) N − i + N (cid:19) . Gini Index is a measure for quantifying sparsity of a distribution. It takes values that rangesfrom 0 to 1. Values closer to 0 is an indication of low sparsity (homogeneity), while values closerto 1 corresponds to higher sparsity (heterogeneity).Thus, we can derive a Gini index for each frequency band in each component of W . We willuse this later, weighted by the relative power of particular IMFs to assess if all channels, or groupsof channels contribute to certain IMFs. For each patient, we quantified the difference in IMF between pairs of seizures for each IMF. Toachieve that, we first computed the product W × IM F i ( t ) with IM F i , representing the multi-dimensional i -th IMF ( k × T matrix). The product yields the matrices X (cid:48) IMF i for all i = 1 , . . . , M timescales. X (cid:48) IMF i reconstructs the i -th IMF in the original space of all channels and frequencybands. For each X (cid:48) IMF i , we computed a distance matrix based on the multivariate Euclideandistance of IMF values for each pair of seizures: D i ( r, c ) = || X (cid:48) IMF i ( t r ) − X (cid:48) IMF i ( t c ) || , where t r and t c are the time epochs of the seizure pair’s onset. Therefore, we obtain M IMF seizure distancematrices per patient, each representing the pairwise seizure distance for a specific IMF.7ote that any seizure-induced changes in the band power will only be present in a few epochs(as we use 30 s epochs). Therefore, the seizures are considered to influence only the fastest IMFs(highest-frequency fluctuations), while they have little effect on the slower IMFs. SupplementaryMaterials additionally shows that our main results reproduce for IMF seizure distances obtainedfrom one epoch before the seizure onset epoch ( t r − and t c − ). In our previous work, we introduced a quantitative measure of how dissimilar two seizures arewithin a patient [2]. Briefly, each epileptic seizure in a patient is analysed in terms of its evolu-tion/pathway through the space of functional network dynamics. The pathways are then comparedto each other using dynamic time warping, accounting for different seizures (or parts of seizures)that are very similar in dynamics but not in duration. The average distance between the warpedseizures is then taken as the dissimilarity measure. For each patient, we can thus obtain a seizuredissimilarity matrix, which captures the pairwise dissimilarity between seizures of that given pa-tient.
In the final part of our analysis, we will investigate if IMF seizure distances can explain seizuredissimilarity in each patient. This part of the analysis was performed for patients with more than5 recorded seizures. We used a linear regression framework, where the seizure dissimilarity is theresponse variable and the IMF seizure distances are the explanatory variables. The observationsare the upper/lower triangle of the matrices representing distances and dissimilarities of uniqueseizure pairs. We also included the temporal distances of seizures (how far apart in time each pairof seizures occurred) as an additional explanatory variable. All variables were standardised beforeincluding them in the model.As a variable selection step for our analysis, we used LASSO [23], which is a sparse shrinkagemethod. Linear regression coefficients were calculated based on least squares, subject to the L penalty. The LASSO also accounted for any collinearity issues between variables. As we wereinterested in detecting positive relationships between the response variable and the explanatoryvariables, we used a constrained positive LASSO, that is, coefficients are constrained to be non-negative. The shrinkage parameter λ for the LASSO model was selected using a 10-fold crossvalidation method from a range of values λ = 10 − , − . . . . , . , (see Suppl. Fig.).After selecting a small number of explanatory variables, an ordinary least squares regressionwas performed, for each patient, to obtain R and confidence intervals for the coefficients. In order to assess if the “best” models selected for each patient could have been found by chance,we performed a permutation test for the adjusted R .In each permutation iteration, we first randomly shifted the seizure onset times within therecording period available. Then, for each iteration, we obtained new IMF seizure distance matricesand we performed the LASSO and linear regression, as described in the previous section, leavingthe response variable unchanged. Finally, we calculated the adjusted R for each iteration. Acrossall iterations, the adjusted R values are used to form an estimation of the distribution of the teststatistic used in the permutation test. P-values were then calculated as the percentage of adjusted8 values that were larger in the permutation distribution. Statistical significance was determinedbased on a significance level of .We also performed an alternative permutation test, where we permuted the order of the seizureswithout permuting the seizure timing. The results are highly comparable and shown in detail inSupplementary Materials. The long-term iEEG recordings for all patients are available at http://ieeg-swez.ethz.ch/ underthe section “Long-term Dataset and Algorithms” [13].Initial signal processing was performed using Matlab version 2019a and Matlab’s built-in func-tions. NMF and MEMD were implemented using the following publicly available functions:• Non-negative matrix factorisation was conducted using the
NNSVD-LRC function from https://sites.google.com/site/nicolasgillis/code [24].• Multivariate EMD was applied using code from [25].For the remainder of the analysis and the construction of all figures we used Python version3.5. Either standard functions obtained from published libraries supported by Python were usedor custom code written in Python. The main functions used in the analysis are listed below:•
Hilbert transform : scipy.signal.hilbert • Mantel test : We used the test function available from http://jwcarr.github.io/MantelTest/ [26].•
LASSO : sklearn.linear_model.Lasso • k-fold cross-validation : sklearn.model_selection.kFold For fitting linear regression models, we used the lm function from stats package implementedin RStudio version 1.1.463 using R version 3.6.3.Our analysis code and data (post processing) can be found on TBC . We analysed the power of multi-channel iEEG signals in the five main frequency bands in termsof their underlying temporal fluctuations on various time-scales for 18 subjects with focal epilepsy.We specifically investigated if fluctuations on particular time scales were driven by particular iEEGfrequency bands, or spatially localised activity. We then explored if these temporal fluctuationscan explain how seizures change in their temporal evolution from one seizure to the next [2].
Following extraction of signal power in the main frequency bands ( δ, θ, α, β, γ ) in 30 s sliding win-dows (no overlap) for each iEEG channel (Fig. 1a,b), we performed dimensionality reduction usingnon-negative matrix factorisation (NMF). NMF effectively groups channels and frequency bandstogether to form components (weights for each component are shown in matrix W , Fig. 1c). The9trength of these components at each time point is then given by the H matrix, which essentiallyyields a time series for each component (Fig. 1d). The original data matrix (Fig. 1b) can beapproximately reconstructed from W ∗ H . The set of time series in H show the fluctuations ofparticular band powers over time.We then determined the different time-scales of fluctuations that were present in each NMFcomponent of H for each patient using Multivariate Empirical Mode Decomposition (MEMD).Fig. 2a shows the MEMD results for a single NMF component in example patient ID06, yielding16 Intrinsic Mode Functions (IMFs) in this case. Faster IMFs (e.g., IMF1, 2 and 3) are oftenthought to contain “noise”, and we investigate this in more detail in Supplementary Materials.However, for simplicity, we will retain all IMFs for the subsequent main results, and refer thereader to Supplementary Materials for a more detailed discussion.Using the instantaneous frequency and amplitude through the Hilbert transform, we can obtaina representation of the marginal spectral densities of each IMF in each dimension/component.Fig. 2b shows the marginal spectral densities summed across all dimensions for each IMF (bluelines) for example patient ID06. Some distinct peaks are seen especially in the slower IMFs, e.g.IMF13 (at ≈ cycle/day), IMF14 (at ≈ . cycle/day), IMF9 ≈ cycles/day), IMF8 ≈ cycles/day), etc. Note that EMD and MEMD as a method essentially acts as a dyadic filter bank[27, 28, 29], thus the dyadic pattern seen in the faster IMFs is not surprising.Across all subjects, fluctuations with a 24h cycle are consistently seen (Fig. 2c). Fluctuationson other time scales can be subject-specific. For 12 out of 18 subjects (ID01, ID02, ID04, ID05,ID07, ID08, ID09, ID11, ID12, ID13, ID17 and ID18) the 24h cycle had the highest density. For5 subjects (ID03, ID06, ID14, ID15 and ID16) the 24h cycle was slightly lower in density, as thehighest density was linked to slower IMFs. For 1 subject (ID10) faster IMFs were marginally higherin density/energy followed by the 24h cycle IMF. Following the observation of a 24h cycle in all subjects, we want to assess the contribution of eachiEEG main frequency band to this 24h fluctuation. To achieve this, we first determined the 24hIMF, which is IMF 13 in example patient ID06 (Fig. 3a). We then calculated the relative powerin each dimension/component of the IMF. For example, in patient ID06, the majority of its poweris concentrated in dimension/component 1 (54%). We also note that the 24h fluctuation doesnot follow the same phase in all dimensions/components of the IMF, potentially indicating thepresence of multiple processes fluctuating on a 24h time scale (see Discussion for more details).For the sake of simplicity, we decided to assess the contribution of different frequency bands overall components/dimensions next.From the dimensionality reduction step, we already know the weight of each iEEG frequencyband for each channel (matrix W, see Fig. 3b). Thus, a summed weight can be obtained for eachfrequency band (across all channels). Finally, a sum weighted by the relative power in the IMFover all dimensions/components can then indicate the relative contribution of each frequency bandto the IMF (Fig. 3c). For most subjects, δ band power contribution is slightly higher compared tothe other frequency bands for the 24h cycle IMF. However, other frequency bands also contributeto the 24h cycle IMF in most patients. We also found that other IMFs have a relatively evencontribution from all frequency bands in all patients (see online data TBC).10 .3 Groups of channels contribute to IMFs slower than 24h Analogously to above, we also investigated the contribution of each channel (in each frequencyband) to an IMF. Specifically, we investigated if the contributions were similar across all channels,or if the contributions were more heterogeneous. We used the weighted Gini index as a measure ofspatial heterogeneity, where 0 (1) indicates a completely homogeneous (heterogeneous) contributionfor each IMF.Fig. 4 shows the distribution of weighted Gini indices in the δ band across all patientsand their IMFs grouped by the IMF dominant timescale of fluctuation. Results for other frequencybands are similar and shown in Supplementary Fig. Overall, the Gini indices are low for all IMFs,indicating that IMFs are not driven by one or few solitary channels. However, there is a cleartendency for slower IMFs (below 1 cycle per day) to display a higher Gini index, indicating thatslower fluctuations may be driven by groups of channels.In summary, we confirmed that the band power IMF fluctuations derived from MEMD areconsistent across patients in terms of a circadian cycle, that the fluctuations are not solely domi-nated by one frequency band (which could indicate strong noise in a particular frequency band),and that the fluctuations are not dominated by one or few channels (another possible indicationof strong noise in one channel). As the final part of our analysis, we investigated if these fluctuations on different time scalesdetected in the IMFs could explain our recent observation of changes in seizure properties overtime in individual patients [2]. Particularly, we previously showed that EEG network dynamicschange from one seizure to the next in every patient, and that these changes could be explainedby hypothetical circadian or longer timescale modulators. Here, we explored if the patient-specificfluctuations detected by the IMFs explain changes in seizure dynamics.For each IMF in each patient, we first determined their corresponding seizure IMF Euclideandistance matrix (Fig.5a,b). For example, in patient’s ID06 IMF 6, we calculated the Euclideandistance of every time point to the time point of the first seizure (Fig.5a) across all dimen-sions/components. By reading out all the Euclidean distances to all the other seizure time points,we obtained the first row of the seizure IMF Euclidean distance matrix (Fig.5b). The same pro-cess can be repeated for all other seizures in this patient. This distance matrix has dimensions ofnumber of seizures by number of seizures and represents how different the IMF state is for eachseizure pair.By using the same techniques as in [2], we obtained a seizure dissimilarity matrix, whichexpresses how dissimilar each pair of seizure pathways are through the space of network dynamics(Fig.5c,d). The seizure dissimilarity matrix represents how each pair of seizures differ withina patient. By relating the seizure dissimilarity for each seizure pair to its corresponding IMFEuclidean distance, we can investigate if there is an association between the two. If the twomeasures are associated, as shown in the example in Fig.5(e), it indicates that the differences inseizure network dynamics can be explained by the IMF state differences.To generalise this approach to all IMFs in a patient, we fitted a multiple linear regressionmodel, where the seizure IMF distances are explanatory variables, and the seizure dissimilarityis the response variable (Fig.6a). We also included the temporal distance between seizures as anexplanatory variable (i.e. how far apart in time each seizure pair occurred). The observations arepairs of seizures, and we fitted a model for each subject. To achieve adequate sample size, we only11ncluded patients with more than 5 seizures. After LASSO variable selection and linear regres-sion, the estimated regression coefficients for subject ID06 are shown in Fig.6(c). The strongestexplanatory effect was seen in the slowest IMF (IMF16) followed by some faster IMFs (IMF3, 4,5 and 6). For this particular subject, according to the model, . of the variability of seizuredissimilarity is explained by explanatory variables (i.e. Adjusted R = 0 . ).Out of all eight patients, six had an Adjusted R around or above 0.6 (Fig.6d). SupplementaryFig. additionally shows that the adjusted R values would have not occurred by chance, exceptfor in ID10. Across all subjects, IMFs around, or slower than, 24h were part of the explanatoryvariables, except ID10 (Fig.6d). However, note that for patient ID10, the model had the lowestadjusted R compared to all the other patients. Temporal distance between seizures remained as anexplanatory variable in three subjects. Faster than 24h IMFs also tended to remain as explanatoryvariables in the models. Overall, a subject-specific combination of different fluctuations appearsto provide a good explanation of seizure variability in most patients.12 omp3IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF8IMF9IMF10IMF11IMF12IMF13IMF14IMF15 IMF16 A m p li t ud e Time (days) (a)
MEMD decomposition of a signal (b) All patients(c) ID06 • All IMFs
ID06 • Comp3
ID02ID01ID03ID04ID05ID08ID07
ID06
ID09ID10ID11ID12ID13ID14ID15ID16ID17ID18
Marginal Hilbert spectrum of fluctuations in total for each patient -3 -2 -1 fluctuations (cycles/day) Marginal Hilbert spectrum of fluctuations of all analytical IMF signals and in total -3 -2 -1 fluctuations (cycles/day) P o w e r IMF2IMF3IMF4IMF5IMF6IMF7IMF8IMF9IMF10IMF11IMF12IMF15IMF13IMF14Total -1 -2 -3 -4 -5 -6 -7 -8 -9 Figure 2: MEMD detects fluctuations on different time-scales across band powerpatterns (components) for each patient. (a) MEMD yields 16 IMFs in example patientID06. Only one dimension of the IMF (corresponding to one NMF component) is shown forsimplicity. IMFs are presented in descending order. (b) Marginal Hilbert spectrum of instantaneousfrequencies for all IMFs (aggregated across all dimensions) in example patient ID06. The blackline represents the Marginal Hilbert frequency spectrum across all IMFs; the last IMF was omittedfor simplicity, as it contains the overall trend. (c) Marginal Hilbert frequency spectrum across all(but last) IMFs for each subject. 13 a) Components
ID06 • 24h Cycle (d)
All patients • 24h Cycle
IMF 13 Time [days] W - component weights (b)
Relative power: D e l t a T he t a A l pha B e t a G a mm a ch1...ch32 I D I D I D I D I D I D I D I D I D I D I D I D I D I D I D I D I D I D D e l t a T he t a A l pha B e t a G a mm a F r e qu e n cy B a nd s PatientsRelative contribution of each iEEG main frequency band to the 24h IMF
Components D e l t a T he t a A l pha B e t a G a mm a (c) Sum of component weights
Figure 3: Contribution of iEEG main frequency bands to the 24h IMF. (a) IMF 13 in ex-ample patient ID06 shows fluctuations on the scale of 24h across all three dimensions/components.Component 1 shows the highest relative power in this IMF. (b) W component weight matrix (sameas Fig. 1c). (c) The sum of the component weights across all channels within each frequency band.(d) Contribution of each iEEG main frequency band to the 24h IMF across all patients obtainedby forming the sum over the matrix in (c) weighted by the relative power in (a). To be ableto compare subjects to each other, each column here has been normalised to form a percentagecontribution. 14 .050.100.150.200.25 δ Frequency Band W e i gh t e d G i n i I nd ex (- . , - . ] (- . , - . ] (- . , - . ] (- . , - . ] (- . , - . ] (- . , . ] ( . , . ] ( . , . ] ( . , . ] ( . , . ] ( . , . ] ( . , . ] All patients
Dominant IMF fluctuation (Cycles/day) - log scale Figure 4: Weighted Gini Index for the δ frequency band across all patients. Acrossall patients, we grouped IMFs based on their dominant IMF fluctuation timescale and show thedistribution of corresponding weighted Gini indices as a violin plot with enclosed box plot. Thewhite dot represents the median, the thick grey bar represents the inter-quartile range.15 e i z u r e d i ss i m il a r i t y Seizures S e i z u r es (a) Euclidean distance from Seizure 1 based on IMF 6 (b)
Seizure dissimilarity (d) (e)
Seizure IMF distanceSeizure IMF distanceSeizures S e i z u r es E u c li d ea n d i s t a n ce ( L ² no r m ) (c) Seizure pathways
ID06 • IMF6 rho = 0.58
Time (days)
Figure 5: Relating seizure dissimilarity and IMF seizure distance.
Throughout this figurewe use example subject ID06 and IMF 6. (a) Euclidean distance of all time points to the firstseizure in terms of IMF 6. Blue dashed vertical lines indicate seizure timing. Dots mark the valueof the IMF distance to the first seizure and colours of dots correspond to the colour scale in (b).(b) Seizure IMF distance matrix for IMF 6. The first row is a representation of the data in (a).(c) Functional network evolution of all seizures projected into 2D space using multidimensionalscaling (MDS) for visualisation of seizure pathways (see [2] for details). Similar seizures tend to beplaced closer together in this projection. Seizures are displayed with distinct colours to distinguishseizure events. The starting points of seizures are marked with a black outline circle. (d) Seizuredissimilarity matrix, capturing the differences in seizure pathways over time between each pairof seizures. (e) Scatter plot of seizure dissimilarity and the seizure IMF distance (Spearman’scorrelation, rho = 0.58, p=0.04). P-value obtained from the Mantel test, and it is FDR correctedfor all IMFs in all patients. 16 istances for all IMFs (a) Response Variable (b) Explanatory Variables Standardised Seizure IMF Distance(c)
All patients (d) Seizure S e i z u r e Seizure SeizureSeizure S e i z u r e StandardisedSeizure Dissimilarity S e i z u r e S e i z u r e −1.2−0.60.00.61.2 Seizure S e i z u r e Standardised Temporal Seizure Distance
ID06ID06 −2.0−1.00.01.02.0ID04 ID06 ID08 ID09 ID10 ID12 ID13 ID140.00.20.40.60.8 A d j u s t e d R ² ID04 ID06 ID08 ID09 ID10 ID12 ID13 ID14
Patients -3 -2 -1 D o m i n a n t F l u c t u a t i on S t a nd a r d i se d r e g r ess i on c o e ff i c i e n t * * * * * ** significance in permutation test C yc l es / d ay * time dist (Intercept)IMF3IMF6IMF5IMF4IMF16−0.25 0.00 0.25 0.50 Standardised regression coefficient E x p l a n a t o r y V a r i a b l es −3.0−2.0−1.00.01.0 Figure 6: IMF seizure distances can explain seizure dissimilarity in most patients. (a) Standardised seizure dissimilarity matrix (response variable). Only the lower triangle of thesymmetric matrix is shown, where each entry serves as an observation. (b) The matrix on theleft shows the standardised temporal seizure distance. Each entry corresponds to the absolutetime difference between seizures. The remaining matrices are standardised seizure IMF distancematrices. (c) Coefficient estimates (black dots), based on ordinary least squares regression forID06, with lines indicating confidence intervals. Only five explanatory variables (IMFs) wereleft after performing variable selection based on constrained LASSO. (d) Summary across patientsbased on OLS models with explanatory variables obtained by the constrained LASSO. Top: Barchart of the Adjusted R . Red stars indicate p − values ≤ . in permutation tests. Bottom: Dotplot indicating the OLS coefficient estimates for the time distance (when this variable remained inthe model) together with OLS coefficient estimates at the corresponding value of dominant IMFfluctuation timescale for each subject. 17 Discussion
We performed a quantitative analysis of multi-channel continuous long-term iEEG recordings from18 focal epilepsy patients. We extracted the band power time series in each channel and identifiedthe temporal fluctuations on various time-scales. Consistent with other studies, we observed thepresence of a circadian rhythm in band power fluctuations across all subjects. Fluctuations on othertime scales were also apparent, but varied between patients, suggesting that various patient-specificfluctuations coexist. In most patients, these patient-specific fluctuations on different timescalesprovided a good explanation for the patient’s seizure dissimilarity (a measure that was previouslyproposed [2] to assess how much the seizure dynamics changed from one seizure to the next withineach patient).Fluctuations on various time-scales of the continuous EEG have been reported in several studiesusing prolonged iEEG recordings. The prevalence of a strong circadian rhythm in EEG patternsis well known since 1969 [30, 31, 32, 33]. Weaker ultradian rhythms have also been reported inlong-term EEG band power during daytime wakefulness in healthy subjects [34, 35]. For example,ultradian modulations are reported to occur with a cycle of 2 and 4 hours in the 9-11 Hz and11-13 Hz frequency bands, respectively [34]. Using a functional connectivity measure of EEG, [8]also reported ultradian fluctuations additional to a strong circadian rhythm. Recently, patient-specific multidien (multi-day) rhythms have been detected in the rate of interictal epileptiformactivity in epilepsy patients [12, 10]. Similarly, daily and longer cycles were shown in the varianceand autocorrelation of EEG signals [36]. Our observations are in agreement with these previousfindings. Apart from the 24h cycle that was apparent in all subjects, we observed fluctuations onultradian and infradian time-scales that were subject-specific.We additionally made two observations about these fluctuations on different time scales. First,we saw that different frequency bands appear to contribute to a similar degree to the circa-dian fluctuation of iEEG band power, although subject-specific patterns of contribution werealso noted. However, our analysis was performed on the multidimensional IMFs (across all di-mensions/components) corresponding to the circadian fluctuation in each patient. The differentcomponents of the IMF were not identical (see e.g. Fig. 3a), indicating that different circadianfluctuations (with different phases) exist in each patient, as has been reported before [33]. Futurework may wish to investigate the different dimensions of the IMF separately in terms of whatfrequency bands contribute to them, and to relate them to other physiological variables such asbody temperature or plasma melatonin [33].The second observation we made was that slower (infradian) fluctuations tended to arise fromgroups of channels, whereas faster (circadian and ultradian) fluctuations tended to arise as a moreequal contribution from all channels. Again, similar to our first observation, we also performed ouranalysis over all dimensions of the IMFs, not distinguishing between IMF components. Individualcomponents of IMFs may well (and do) have spatially heterogeneous contributors. To fully uncoverthe spatial and frequency band contributions to each component of each IMF, we suggest thatfuture work performs an iterative dimensionality reduction and empirical mode decomposition tooptimally find components and component contributions for each IMF.Our main goal was to investigate if fluctuations in long-term iEEG band power explain seizuredissimilarity within each patient. Seizure dissimilarity is a measure that was proposed to assesshow different any pair of seizures are in a given patient in terms of their seizure evolution anddynamics. Previous work has also shown that seizure dissimilarity was well-explained by processesincorporating Gaussian noise, circadian and/or slower timescales changes in most patients [2]. Inagreement, we found that circadian or infradian fluctuations contributed strongly in most patients18n explaining seizure dissimilarity. Interestingly, we also found many faster (ultradian) fluctuationsas explanatory variables in most patients. These could be contributing explanatory power throughwhat previously was modelled as noise. With larger datasets using many more seizures recordedover a longer period, future work should investigate ultradian contributions carefully and assess ifnoise would perform as well as the cumulative ultradian contributions.We also wish to point out that the band power fluctuations did not account for all the variabilityof seizure dissimilarity in most patients. The highest R was around 0.8 and the unexplainedvariability based on the models suggests that there are additional factors that impact seizuredynamics. One such factor may be the anti-epileptic drug (AED) level at any given time. Patientsundergoing such prolonged iEEG monitoring often undergo AED withdrawal as part of the clinicalevaluation. It is also well-known that AED changes and withdrawal can change the severity anddynamics of seizures. For example, bilateral tonic-clonic seizures are more prevalent when AEDlevels are reduced [37]. We were unable to include this information in the current study, but futurestudies may wish to investigate AED levels as additional potential explanatory variables for seizuredissimilarity.In a more general context, our work is another contribution to the wider literature of explainingictal features from interictal EEG features, or hypothesised circadian/infradian rhythms. Mostnotably, the ictal feature of interest is often the seizure timing (when a seizure occurs). For example,studies have established that there is often a (subject-specific) relationship between fluctuationsof interictal EEG features and the timing of ictal events [12, 10, 8]. A number of studies alsodemonstrated that seizure occurrence is underpinned by a circadian/infradian rhythms [12, 10, 38,39, 31]. Interestingly, we found no evidence of an association between band power fluctuationsof the interictal EEG and seizure occurrence (data not shown). Seizures were not statisticallymore likely to occur during particular phases of particular IMFs in most patients in our dataset. This is in agreement with a previous study [8] that reported functional network fluctuations,rather than band power fluctuations, to be more predictive of seizure timing. Therefore, futurework should investigate temporal fluctuations in a range of EEG features, such as band power,functional connectivity [8], high frequency oscillations [9], variance and autocorrelation [36], justto name a few. Apart from seizure timing, our work has shown that these fluctuations on differenttimescales of the interictal EEG also explain other features of the seizures (in our case seizuredissimilarity). Future work should explore this avenue further to elucidate the exact processes andtimescales that modulate or dictate the various aspects of a seizure.A variety of possible factors, both external and internal, could explain the observations aboveof interictal EEG fluctuations on different time scales, and correlated changes in seizure features.External factors may be environmental processes such as lower atmospheric pressure, high airhumidity and temperature [40] or processes related to nutrition, such as caffeine [41].Internal factors, linked to external factors, include emotional/physical stress [42], low qualityof sleep [43], and fatigue [44]. Other internal factors include hormonal rhythms (e.g. menstrualcycles) [45, 46] and circadian rhythms [47]. Seizure frequency, seizure dynamics, interictal epilepticphenomena (e.g. spikes, bursts, HFOs) and general interictal EEG patterns all might be affectedby these time-varying factors [48, 49]. However, the exact relationship between modulators andepileptic processes remains to be explored in future.In conclusion, seizure occurrence/timing may not be the only aspect of a seizure that is corre-lated with fluctuating interictal EEG features. Our work has shown that the EEG evolution anddynamics during a seizure may be modulated by such interictal fluctuations on infradian, circadian,and ultradian time-scales. The implication is that EEG fluctuations on different timescales maypredict seizure dynamics. Due to the tremendous impact epilepsy has on the quality of a patient’s19ife, it is necessary to provide predictions of various seizure features based on the patient’s needs,including seizure severity. Furthermore, if interictal fluctuations can be manipulated, then theymay also causally alter seizure dynamics. As a novel therapeutic approach, seizures may then berendered less severe as a consequence. References [1] C. E. Stafstrom and L. Carmant, “Seizures and Epilepsy: An Overview for Neuroscientists,”
Cold Spring Harb Perspect Med , vol. 5, p. a022426, June 2015.[2] G. M. Schroeder, B. Diehl, F. A. Chowdhury, J. S. Duncan, J. d. Tisi, A. J. Trevelyan,R. Forsyth, A. Jackson, P. N. Taylor, and Y. Wang, “Seizure pathways change on circadianand slower timescales in individual patients with focal epilepsy,”
PNAS , May 2020.[3] V. K. Jirsa, W. C. Stacey, P. P. Quilichini, A. I. Ivanov, and C. Bernard, “On the nature ofseizure dynamics,”
Brain , vol. 137, pp. 2210–2230, Aug. 2014.[4] M. L. Saggio, D. Crisp, J. Scott, P. J. Karoly, L. Kuhlmann, M. Nakatani, T. Murai, M. Düm-pelmann, A. Schulze-Bonhage, A. Ikeda, M. Cook, S. V. Gliske, J. Lin, C. Bernard, V. Jirsa,and W. Stacey, “Epidynamics characterize and navigate the map of seizure dynamics,” bioRxiv ,p. 2020.02.08.940072, Feb. 2020.[5] B. Oken and K. Chiappa, “Short-term variability in EEG frequency analysis,”
Electroen-cephalography and Clinical Neurophysiology , vol. 69, pp. 191–198, Mar. 1988.[6] C. Geier, K. Lehnertz, and S. Bialonski, “Time-dependent degree-degree correlations in epilep-tic brain networks: from assortative to dissortative mixing,”
Front. Hum. Neurosci. , vol. 9,2015.[7] C. Geier and K. Lehnertz, “Long-term variability of importance of brain regions in evolvingepileptic brain networks,”
Chaos , vol. 27, p. 043112, Apr. 2017.[8] G. D. Mitsis, M. N. Anastasiadou, M. Christodoulakis, E. S. Papathanasiou, S. S. Papacostas,and A. Hadjipapas, “Functional brain networks of patients with epilepsy exhibit pronouncedmultiscale periodicities, which correlate with seizure onset,”
Human Brain Mapping , vol. 41,no. 8, pp. 2059–2076, 2020.[9] S. V. Gliske, Z. T. Irwin, C. Chestek, G. L. Hegeman, B. Brinkmann, O. Sagher, H. J. L.Garton, G. A. Worrell, and W. C. Stacey, “Variability in the location of high frequencyoscillations during prolonged intracranial EEG recordings,”
Nat Commun , vol. 9, p. 2155,Dec. 2018.[10] P. J. Karoly, D. R. Freestone, R. Boston, D. B. Grayden, D. Himes, K. Leyde, U. Senevi-ratne, S. Berkovic, T. O’Brien, and M. J. Cook, “Interictal spikes and epileptic seizures: theirrelationship and underlying rhythmicity,”
Brain , vol. 139, pp. 1066–1078, Apr. 2016.[11] E. C. Conrad, S. B. Tomlinson, J. N. Wong, K. F. Oechsel, R. T. Shinohara, B. Litt, K. A.Davis, and E. D. Marsh, “Spatial distribution of interictal spikes fluctuates over time andlocalizes seizure onset,”
Brain , vol. 143, pp. 554–569, Feb. 2020.2012] M. O. Baud, J. K. Kleen, E. A. Mirro, J. C. Andrechak, D. King-Stephens, E. F. Chang, andV. R. Rao, “Multi-day rhythms modulate seizure risk in epilepsy,”
Nature Communications ,vol. 9, p. 88, Jan. 2018.[13] A. Burrello, L. Cavigelli, K. Schindler, L. Benini, and A. Rahimi, “Laelaps: An Energy-Efficient Seizure Detection Algorithm from Long-term Human iEEG Recordings without FalseAlarms,” in ,(Florence, Italy), pp. 752–757, IEEE, Mar. 2019.[14] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factoriza-tion,”
Nature , vol. 401, pp. 788–791, Oct. 1999.[15] S. M. Atif, S. Qazi, and N. Gillis, “Improved SVD-based initialization for nonnegative matrixfactorization using low-rank correction,”
Pattern Recognition Letters , vol. 122, pp. 53–59, May2019.[16] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung,and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinearand non-stationary time series analysis,”
Proc. R. Soc. Lond. A , vol. 454, pp. 903–995, Mar.1998.[17] N. E. Huang, M.-L. C. Wu, S. R. Long, S. S. P. Shen, W. Qu, P. Gloersen, and K. L. Fan,“A Confidence Limit for the Empirical Mode Decomposition and Hilbert Spectral Analysis,”
Proceedings: Mathematical, Physical and Engineering Sciences , vol. 459, no. 2037, pp. 2317–2345, 2003.[18] A. Y. Kaplan, A. A. Fingelkurts, A. A. Fingelkurts, S. V. Borisov, and B. S. Darkhovsky,“Nonstationary nature of the brain activity as revealed by EEG/MEG: Methodological, prac-tical and conceptual challenges,”
Signal Processing , vol. 85, pp. 2190–2212, Nov. 2005.[19] A. A. Fingelkurts and A. A. Fingelkurts, “Operational Architectonics of the Human BrainBiopotential Field: Towards Solving the Mind-Brain Problem,”
Brain and Mind , vol. 2,pp. 261–296, Dec. 2001.[20] G. Rilling, P. Flandrin, P. Goncalves, and J. M. Lilly, “Bivariate Empirical Mode Decompo-sition,”
IEEE Signal Processing Letters , vol. 14, pp. 936–939, Dec. 2007.[21] Y. Lv, R. Yuan, and G. Song, “Multivariate empirical mode decomposition and its applicationto fault diagnosis of rolling bearing,”
Mechanical Systems and Signal Processing , vol. 81,pp. 219–234, Dec. 2016.[22] N. P. Hurley and S. T. Rickard, “Comparing Measures of Sparsity,” arXiv:0811.4706 [cs,math] , Apr. 2009. arXiv: 0811.4706.[23] R. Tibshirani, “Regression Shrinkage and Selection via the Lasso,”
Journal of the Royal Sta-tistical Society. Series B (Methodological) , vol. 58, no. 1, pp. 267–288, 1996.[24] A. M. Syed, S. Qazi, and N. Gillis, “Improved SVD-based Initialization for Nonnegative MatrixFactorization using Low-Rank Correction,”
Pattern Recognition Letters , vol. 122, pp. 53–59,May 2019. arXiv: 1807.04020. 2125] N. Rehman and D. P. Mandic, “Multivariate empirical mode decomposition,”
Proc. R. Soc.A. , vol. 466, pp. 1291–1302, May 2010.[26] N. Mantel, “The Detection of Disease Clustering and a Generalized Regression Approach,”
Cancer Res , vol. 27, pp. 209–220, Feb. 1967.[27] Z. Wu and N. E. Huang, “A study of the characteristics of white noise using the empiricalmode decomposition method,”
Proc. R. Soc. Lond. A , vol. 460, pp. 1597–1611, June 2004.[28] P. Flandrin, G. Rilling, and P. Goncalves, “Empirical mode decomposition as a filter bank,”
IEEE Signal Processing Letters , vol. 11, pp. 112–114, Feb. 2004.[29] N. ur Rehman and D. P. Mandic, “Filter Bank Property of Multivariate Empirical ModeDecomposition,”
IEEE Trans. Signal Process. , vol. 59, pp. 2421–2426, May 2011.[30] H. Scheich, “Interval histograms and periodic diurnal changes of human alpha rhythms,”
Electroencephalogr Clin Neurophysiol , vol. 26, p. 442, Apr. 1969.[31] D. C. Spencer, F. T. Sun, S. N. Brown, B. C. Jobst, N. B. Fountain, V. S. S. Wong, E. A.Mirro, and M. Quigg, “Circadian and ultradian patterns of epileptiform discharges differby seizure-onset location during long-term ambulatory intracranial monitoring,”
Epilepsia ,vol. 57, pp. 1495–1502, Sept. 2016.[32] M. K. Smyk and G. van Luijtelaar, “Circadian Rhythms and Epilepsy: A Suitable Case forAbsence Epilepsy,”
Front. Neurol. , vol. 11, 2020.[33] D. Aeschbach, J. R. Matthews, T. T. Postolache, M. A. Jackson, H. A. Giesen, and T. A.Wehr, “Two circadian rhythms in the human electroencephalogram during wakefulness,”
American Journal of Physiology-Regulatory, Integrative and Comparative Physiology , vol. 277,pp. R1771–R1779, Dec. 1999.[34] D. A. Kaiser, “Ultradian and Circadian Effects in Electroencephalography Activity,” p. 4,2008.[35] F. Chapotot, C. Jouny, A. Muzet, A. Buguet, and G. Brandenberger, “High frequency wak-ing EEG: reflection of a slow ultradian rhythm in daytime arousal,”
NeuroReport , vol. 11,pp. 2223–2227, July 2000.[36] M. I. Maturana, C. Meisel, K. Dell, P. J. Karoly, W. D’Souza, D. B. Grayden, A. N. Burkitt,P. Jiruska, J. Kudlacek, J. Hlinka, M. J. Cook, L. Kuhlmann, and D. R. Freestone, “Criticalslowing down as a biomarker for seizure susceptibility,”
Nat Commun , vol. 11, May 2020.[37] M. C. Pensel, M. Schnuerch, C. E. Elger, and R. Surges, “Predictors of focal to bilateral tonic-clonic seizures during long-term video-EEG monitoring,”
Epilepsia , vol. 61, no. 3, pp. 489–497,2020.[38] P. J. Karoly, H. Ung, D. B. Grayden, L. Kuhlmann, K. Leyde, M. J. Cook, and D. R. Freestone,“The circadian profile of epilepsy improves seizure forecasting,”
Brain , vol. 140, pp. 2169–2182,Aug. 2017. 2239] P. J. Karoly, D. M. Goldenholz, D. R. Freestone, R. E. Moss, D. B. Grayden, W. H. Theodore,and M. J. Cook, “Circadian and circaseptan rhythms in human epilepsy: a retrospective cohortstudy,”
The Lancet Neurology , vol. 17, pp. 977–985, Nov. 2018.[40] F. Rakers, M. Walther, R. Schiffner, S. Rupprecht, M. Rasche, M. Kockler, O. W. Witte,P. Schlattmann, and M. Schwab, “Weather as a risk factor for epileptic seizures: A case-crossover study,”
Epilepsia , vol. 58, no. 7, pp. 1287–1295, 2017.[41] R. R. van Koert, P. R. Bauer, I. Schuitema, J. W. Sander, and G. H. Visser, “Caffeine andseizures: A systematic review and quantitative analysis,”
Epilepsy Behav , vol. 80, pp. 37–47,2018.[42] E. Baldin, W. A. Hauser, A. Pack, and D. C. Hesdorffer, “Stress is associated with an increasedrisk of recurrent seizures in adults,”
Epilepsia , vol. 58, no. 6, pp. 1037–1046, 2017.[43] C. W. Bazil and T. S. Walczak, “Effects of Sleep and Sleep Stage on Epileptic and NonepilepticSeizures,”
Epilepsia , vol. 38, no. 1, pp. 56–62, 1997.[44] O.-Y. Kwon, H. S. Ahn, and H. J. Kim, “Fatigue in epilepsy: A systematic review and meta-analysis,”
Seizure , vol. 45, pp. 151–159, Feb. 2017.[45] J. Bauer, “Interactions Between Hormones and Epilepsy in Female Patients,”
Epilepsia , vol. 42,no. s3, pp. 20–22, 2001.[46] E. Taubøll, L. Sveberg, and S. Svalheim, “Interactions between hormones and epilepsy,”
Seizure , vol. 28, pp. 3–11, May 2015.[47] S. Khan, L. Nobili, R. Khatami, T. Loddenkemper, C. Cajochen, D.-J. Dijk, and S. H.Eriksson, “Circadian rhythm and epilepsy,”
The Lancet Neurology , vol. 17, pp. 1098–1108,Dec. 2018.[48] D. R. Freestone, P. J. Karoly, and M. J. Cook, “A forward-looking review of seizure predic-tion:,”
Current Opinion in Neurology , vol. 30, pp. 167–173, Apr. 2017.[49] R. A. B. Badawy, D. R. Freestone, A. Lai, and M. J. Cook, “Epilepsy: Ever-changing statesof cortical excitability,”