Circadian distribution of epileptiform discharges in epilepsy: Candidate mechanisms of variability

Epilepsy is a serious neurological disorder characterised by a tendency to have recurrent, spontaneous, seizures. Classically, seizures are assumed to occur at random. However, recent research has uncovered underlying rhythms both in seizures and in key signatures of epilepsy—so-called interictal epileptiform activity—with timescales that vary from hours and days through to months. Understanding the physiological mechanisms that determine these rhythmic patterns of epileptiform discharges remains an open question. Many people with epilepsy identify precipitants of their seizures, the most common of which include stress, sleep deprivation and fatigue. To quantify the impact of these physiological factors, we analysed 24-hour EEG recordings from a cohort of 107 people with idiopathic generalized epilepsy. We found two subgroups with distinct distributions of epileptiform discharges: one with highest incidence during sleep and the other during day-time. We interrogated these data using a mathematical model that describes the transitions between background and epileptiform activity in large-scale brain networks. This model was extended to include a time-dependent forcing term, where the excitability of nodes within the network could be modulated by other factors. We calibrated this forcing term using independently-collected human cortisol (the primary stress-responsive hormone characterised by circadian and ultradian patterns of secretion) data and sleep-staged EEG from healthy human participants. We found that either the dynamics of cortisol or sleep stage transition, or a combination of both, could explain most of the observed distributions of epileptiform discharges. Our findings provide conceptual evidence for the existence of underlying physiological drivers of rhythms of epileptiform discharges. These findings should motivate future research to explore these mechanisms in carefully designed experiments using animal models or people with epilepsy.

Potential subgroups are identified using k-means clustering, which is executed using the built-in MATLAB function kmeans.In our computation, we repeated the clustering 10 times using new initial cluster centroid positions.
We repeat the clustering using different numbers of clusters (from 1 to 10), and we employ the Calinski-Harabasz criterion (using the MATLAB function evalclusters) to identify the optimal number.This method is based on the Calinski-Harabasz value, which measures the goodness of fit of the data partition: the higher the value, the better that specific number of clusters succeeds in partitioning the data.Moreover, to test the robustness of our choice of considering a 1-hour time window in our study, we modified the bin width of the histograms by defining the total number of bins (N bin ) in constructing these histograms.Setting N bin = 24 corresponds to assessing the distribution of the EDs on an hourly basis as presented in the Manuscript.

Impact of histogram bin width
In the Manuscript, we organized IGE subjects into two clusters based on the similarities in their cumulative epileptiform discharge patterns on an hourly basis.We test the robustness of our result by repeating the clustering for a range of values of N bin (see Table A and Fig B).The size of these two clusters across the different bin width values is similar, which shows the findings in the main manuscript are robust with respect to the specific choice of bin width.

Groups vs type of epilepsy
Table B provides a detailed breakdown of the different syndromes within the two groups.To assess the dependence of the groups on the specific syndrome, we fit a linear model using the MATLAB function fitlm and the model formula Group ∼ Syndrome.
The p-value for the test is p = 0.756, suggesting there is no strong evidence of an association between the variability in EDs occurrences in the two groups and specific syndromes.

Temporal patterns of ED events
To test the dependence of EDs occurrence on the circadian hour time and sleep in Group 1 and 2, we employed a mixed-effect Poisson regression model with the formula: where observed variable ED = (ED 1 , ..., ED N ) is a 107×24 matrix where each ED i corresponds to the ED occurrence for the i th individual of the IGE cohort during the 24-hour time window.The predictor Time represents the circadian time (hours), while Sleep is a 24-element vector whose entrance is 1 if the individual is sleeping and 0 otherwise.Due to the intra-subject variability, we introduce the variable Subject as a random factor.We used the MATLAB function fitlme to estimate the slope of the traces corresponding to each measurement.
Table C summarizes the results of the linear regression for Groups 1 and 2. In both groups, there is a statistically significant change in ED counts across time blocks and sleep (p-value <0.001).This suggests an impact of sleep on the ED occurrence in Group 2 as well as in Group 1.This result can be explained by noticing that the morning peak in ED events recorded in Group 2 starts during sleep time.

Null model vs external perturbation
To test whether the model outputs are only due to random noise, we defined the 'null model' as the model with no external input, i.e., λ ext (t) = 0. We run this model for both groups and we perform a statistical analysis to check whether there is a statistically significant difference across the day between the null model (Null) and the simulation with the external factor (Sim).The linear model was fit using the MATLAB function fitlm with the formula: where Y is either the null model (no external perturbation) or the model simulation with sleep and CORT only for Group 1 and 2, respectively.Tables D and E summarize the results for the null model and the model accounting for the external factors, respectively.For both groups, Sim shows a statistically relevant change during the day (p-time<0.001),while Null shows no changes (p-time>0.1).This result suggests that there is an effect of the external factors on the model simulation and it is not due only to random noise.

R 2 statistic for the combined mechanism
To assess the goodness-of-fit of the model, we could also have considered R 2 , which corresponds to the proportion of variance explained by the model.In order to assess the robustness of the parameter estimations of p S and p C , we considered two other approaches: the maximum-likelihood estimation (MLE) and the Monte Carlo Markov Chain (MCMC).

Maximum-likelihood estimation
The MLE is a method of estimating the parameters of a statistical model given a set of observations by maximising the likelihood function, L(X|θ) (which is effectively equivalent to a cost-function approach), with X as the vector of observed data and θ the parameter vector that we want to estimate.In our study, X is the ED rate and θ = (p S , p C ).We define the vector Y θ as the ED rate estimated by our model for a certain value of θ.We then assume that the ED rate X comes from a normal distribution with mean Y θ and 0.2 standard deviations, that is X ∼ N (Y θ , 0.2).The value 0.2 corresponds to the mean of the standard deviation values computed hourly across the cohort.
We therefore compute the values of L(X|p S , p C ) over a grid where 0 ≤ p S ≤ 1.5 and 0 ≤ p C ≤ 1.5.Fig G illustrates the likelihood values for Group 1 (panel A) and Group 2 (panel B).We find that the likelihood function is maximal when p S = 1.1 and p C = 0.6 for Group 1 and p S = 0 and p C = 1 for Group 2. This result is consistent with what has been found using the RSS approach in the main text.

Monte Carlo Markov Chain
An alternative method to estimate p S and p C is the MCMC.In our study, we implemented an adaptive Metropolis (AM) algorithm to estimate the parameters θ = (p S , p C ) as described in [1,2].
In short, at every sampler iteration i, we suggest a proposal distribution π := N tr (θ i−1 , Σ i ), where θ i−1 is the latest accepted value and Σ i is the covariance matrix defined in Eq. 1, and we extract the proposed parameter θ i from π.We then compute the probability α := min 1, )), the parameter is accepted and the mean and the covariance matrix of the proposal distribution are updated, that is θ i+1 ∼ N tr (θ i , Σ i+1 ).Otherwise, the proposed parameter is rejected and θ i+1 ∼ N tr (θ i−1 , Σ i+1 ).The result is a time series of accepted proposed parameters {θ} i (posterior distribution).
Fig H shows the time series for four different initial values (sampled from the prior distribution π) of the parameters p S and p C for both Group 1 and Group 2 (top row) and the corresponding auto-correlation values (bottom row).We ran the MCMC for 5 thousand iterations of which the first one thousand were discarded as burn-in.From the remaining 4 thousand, we recorded samples every 4 iterations for a total of one thousand.
The final values are given by the median of {θ} i and are reported in Table F.Although the single values are not identical given the different approaches (Bayesian vs Frequentist), the parameters estimated with the RSS and MLE approaches are within the 95% confidence interval estimated with MCMC.More precisely, sleep is shown to be the main driver in Group 1 (p S > 1), although there is a non-negligible CORT component (p C ∼ 0.3).Conversely, in Group 2, sleep is not playing a role in shaping the ED distribution (p C ∼ 0) while CORT is the main driver (p C > 1).

Representative model simulation
The mathematical model used in this study is based on the normal form of a subcritical Hopf bifurcation: żi where z i (t) is dynamics of the i th node (with i = 1, ..., N ), W (t) is a complex Wiener process, λ(t) is the excitability of node i, λ base the baseline level of excitability, λ ext (t) the external perturbations to the excitability, ω determines the oscillation frequency and α determines the level of noise for the complex Wiener process W (t). Typical parameter values for the model are given in Table 2 in the main text, whilst A is an adjacency matrix, i.e.A i,j is 1 if there is a connection between the i th and j th regions and 0 otherwise.For simplicity, all simulations were performed with a directed and connected 4-node graph (N = 4) (Fig 9 in the main text).Numerical simulations were obtained using an Euler-Maruyama scheme to find approximate solutions to the system of stochastic differential equations (SDEs) with dt = 10 −3 .The method was implemented in MATLAB R2021a (MathWorks Inc., Natick, MA) to simulate 24-hour brain activity.The resulting ED rate histogram is reported in (D).
Fig A. Selecting optimal number of clusters within the IGE cohort.Calinski-Harabasz values corresponding to different numbers of clusters and time windows.
We can cluster the patients based on the cross-correlation computed when the time series are aligned to the patient's sleep onset and offset.Fig C illustrates the output when the series are aligned to clock time (column A), sleep onset (column B), and sleep offset (column C).Although some patients classified into Group 1 when the clock time is used are now clustered within Group 2 when sleep times are considered (and vice-versa), the overall distributions are robust.Therefore, we decided to use the clock time in our study (see main text).

Fig
Fig C. Clustering when the time series are aligned to sleep times.ED hourly rate patterns were aligned to the clock time (A), sleep onset (B) and sleep offset (C).A pairwise cross-correlation matrix (107 × 107) was calculated to establish similarities within the IGE cohort (first row).The corresponding Group 1 and 2 are shown in the middle and bottom rows, respectively.

Fig
Fig D. Network selection.Examples of simulations with different network structures (left column) and the corresponding model simulations for Group 1 (middle column) and Group 2 (right column).
Fig F. R 2 values for the combined mechanisms.Values of the R 2 statistic computed over a grid of values of p S and p C for Group 1 (A) and Group 2 (B).
10. Robustness of estimating p S and p C .

Fig G .
Fig G. Likelihood function values for the combined mechanism.Values of the likelihood function (L(X|p S , p C )) computed over a grid at intervals of 0.1 of values of p S and p C for Group 1 (A) and Group 2 (B).

Fig
Fig H. Results for the MCMC algorithm.Time series of the parameter values obtained with the adaptive Metropolis for the two parameters p S and p C for both Group 1 and Group 2 (top row) and the corresponding autocorrelation function (bottom row).The different colours correspond to the different initialization of the parameters.
Fig I. Selecting optimal number of clusters within the two groups.Calinski-Harabasz values corresponding to Group 1 (blue) and Group (red).

Fig J illustrates the
Fig J illustrates the subsequent subgroups.Interestingly, this analysis identifies the subjects in Group 2 responsible for the evening peak, which is unlikely to be explained by the CORT dynamics.

Fig K illustrates the
Fig K illustrates the optimal fit with the corresponding value of p S , p C , and R 2 .

D
Fig M. Simulation representative.Numerical simulations of the variable z i (A) and λ i (B), with i = 1, ..., 4. The λ ext used in the simulation is given by the sum of λ ext,sleep and λ ext,CORT shown in (C).The resulting ED rate histogram is reported in (D).

Table A .
Number of subjects in each group for the different number of bins.
Fig B. Model results compared with IGE data for different numbers of bins.Left column: Histogram of epileptiform discharges from Group 1 with IGE (blue) and histogram of epileptiform discharges simulated using the model with λ ext defined to mimic the different brain excitability during sleep stages (green).Right column: Histogram of epileptiform discharges from Group 2 with IGE (red) and histogram of epileptiform discharges simulated using the model with λ ext defined

Table C .
Summary of the linear mixed-effects modelling of ED occurrence.Summary of the linear mixed-effects modelling of ED occurrence in Group 1. Estimates, 95% confidence intervals, and p-value (bold, p-value <0.05) for each predictor of ED are shown.

Table D .
Summary of the linear modelling of the null model.Summary of the linear mixed-effects modelling of ED occurrence in Group 1 and 2. Estimates, 95% confidence intervals, and p-value (bold, p-value <0.05) for each predictor of ED are shown.

Table E .
Summary of the linear modelling of with perturbation.Summary of the linear modelling of ED occurrence simulated by our model when Sleep is the only external factor in Group 1 and CORT in Group 2. Estimates, 95% confidence intervals, and p-value (bold, p-value <0.05) for each predictor of ED are shown.
B Fig E. Model results compared with IGE data.(A) Histogram of EDs from Group 1 with IGE (blue) and histogram of EDs simulated using the model with λ ext defined to mimic the impact of CORT on the brain excitability (green).(B) Histogram of EDs from Group 2 with IGE (red) and histogram of EDs simulated using the model with λ ext defined to mimic the different brain excitability during sleep stages (green).

Table F .
Parameter estimation with MCMC.Summary of the estimations for the parameters p S and p C obtained with the AM MCMC for each initial condition and each of the two groups, and the corresponding 95% confidence interval identified by the first and third quartiles (2.5% percentile and 97.5% percentile, respectively) of the posterior distribution.