Temporal Evolution of Multiday, Epileptic Functional Networks Prior to Seizure Occurrence

Epilepsy is one of the most common neurological disorders, characterized by the occurrence of repeated seizures. Given that epilepsy is considered a network disorder, tools derived from network neuroscience may confer the valuable ability to quantify the properties of epileptic brain networks. In this study, we use well-established brain network metrics (i.e., mean strength, variance of strength, eigenvector centrality, betweenness centrality) to characterize the temporal evolution of epileptic functional networks over several days prior to seizure occurrence. We infer the networks using long-term electroencephalographic recordings from 12 people with epilepsy. We found that brain network metrics are variable across days and show a circadian periodicity. In addition, we found that in 9 out of 12 patients the distribution of the variance of strength in the day (or even two last days) prior to seizure occurrence is significantly different compared to the corresponding distributions on all previous days. Our results suggest that brain network metrics computed fromelectroencephalographic recordings could potentially be used to characterize brain network changes that occur prior to seizures, and ultimately contribute to seizure warning systems.


Introduction
The human brain is a dynamic system that undergoes several dynamic changes due to, for example, different cognitive processes, states of vigilance or motor tasks [1]. The brains of people with epilepsy (PWE) have an additional dynamic change in their dynamic repertoire, making transitions between "normal activity" and "seizure activity" [2]. The occurrence of repeated seizures is the main characteristic of epilepsy and seizures that cannot be controlled by medications cause a significant burden in the lives of PWE [3,4]. Epilepsy is considered to be a network disorder, with several studies suggesting that even focal seizures arise from large-scale networks [5][6][7][8]. A network can be characterized as being composed of nodes that are connected with edges. In a large-scale brain network, the nodes represent brain regions whilst the edges denote statistical interactions between the nodes (functional network) or anatomical connections between the nodes (structural network) [9]. Hence, to better understand epilepsy it may be valuable to study the topological properties of the epileptic brain networks.
Brain networks can be studied at different spatial and temporal scales using different data modalities [10]. Several studies that aimed to investigate the temporal evolution of epileptic brain networks from minutes to hours and days used intracranial electroencephalographic recordings (iEEG) and built functional networks [11][12][13][14][15]. They assumed each iEEG channel constituted one node and inferred statistical connections between the nodes using linear and non-linear connectivity measures. To characterize topological properties of the brain networks they used metrics from graph theory such as degree (number of links of nodes), betweenness centrality (quantifying the influence that a node has over the flow of information in the network), clustering coefficient (quantifying the tendency of nodes to cluster together) and shortest path length [16,17]. These studies reported high temporal variability of the graph theory metrics across many days with the presence of daily rhythms having a key role in this temporal variation. Moreover, they reported that the brain network evolves from a random structure to a more regular structure during seizures with increasing randomness after seizure termination.
Although iEEG recordings have excellent data quality, their acquisition is highly invasive as it requires brain surgery. In addition, the iEEG electrodes are implanted in spatially limited brain regions and hence do not record electrical activity from the whole brain. In this study, we use longitudinal non-invasive scalp electroencephalographic recordings (EEG) from PWE and build evolving functional networks that span from 2 to 11 days. Using metrics from graph theory we quantify the topological properties of the brain networks and explore their temporal evolution across many days prior to seizure occurrence. We further investigate the presence of periodic patterns in this temporal evolution. Considering the influence of daily rhythms in the temporal evolution of the graph theory metrics, we ask whether the graph metrics that are computed from brain networks of the day prior to seizure occurrence show particular changes in the brain network topology compared to the previous days. Since the EEG acquisition was accompanied with an electrocardiogram (ECG) acquisition, we also investigated the temporal evolution of cardiac metrics (heart rate and heart rate variability). Finally, we discuss the study findings and highlight their importance.

Participants Recruitment
Study participants were recruited from July 2017 to February 2020 via the European RADAR-CNS (Remote Assessment of Disease and Relapse-Central Nervous System) study. During this period, 71 participants with an epilepsy diagnosis were admitted to the Epilepsy Monitoring Unit (EMU) of King's College Hospital, London as part of their clinical workup. The length of the patients' stay in the EMU was no more than two weeks and was determined by clinical considerations. All patients were monitored continuously via a video-EEG system that recorded simultaneously brain and cardiac activity resulting in 21-channel EEG (Fp1, Fp2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, A1, A2, P3, Pz, P4, T6, O1, O2; electrode placement according to the modified Maudsley system) and two ECG signals. The sampling frequency was either 256 Hz or 512 Hz. An expert neurologist (E.B.) reviewed all EEG recordings and annotated the timestamps of the seizures onset and offset. More details regarding the data acquisition can be found in Bruno et al. [18].

Inclusion and Exlusion Criteria
Since we were interested in the multiday temporal evolution of brain networks prior to seizure occurrence we excluded from the analysis all patients who did not experience seizures during their stay in the EMU. In addition, patients who had seizures for which there was not at least a two-day seizure-free period prior to the seizure occurrence were also excluded from this analysis. Hence, the final dataset that we analyzed comprised 12 participants (one seizure per patient) for whom the seizure-free period prior to the seizure occurrence spanned from 2 to 15 days (Table 1).

Ethics Statement
The study was conducted in accordance with the Declaration of Helsinki. All participants gave written informed consent and the study procedures were approved by the London Fulham Research Ethics Committee (16/LO/2209; IRAS project ID216316).

EEG Measurements
We visually inspected the available EEG recordings and selected a 20 s artifact-free epoch from every hour (i.e., 24 EEG epochs/day). Across all patients, the minimum amount of available daily EEG epochs was 21 (missing epochs originated due to the disconnection of the electrodes or low signal quality).

Functional Network
We built functional networks for each 20 sec EEG epoch by associating the network nodes to electrodes, whilst the edges between the nodes and were derived by the phase locking value (PLV) [19][20][21], where is the number of samples and , ( ) is the instantaneous phase difference between the EEG signals that correspond to nodes and . The phases were extracted from the filtered signals using the Hilbert transform. To account for the possible effects of volume conduction [22], connections at zero time lag , = (∑ , ( ) ) were discarded. Further, to account for the presence of connections that are due to chance or due to the finite length of the signals we additionally used surrogate signals. In particular, for every pair of signals, we generated 99 surrogate pairs of signals using the iterative amplitude-adjusted Fourier transform (IAAFT) with 10 iterations [23,24]. Hence, for every PLV value of the ( , ) pair of nodes we obtained a distribution of 99 PLV values from the surrogate signals. If the PLV value of the ( , ) pair was higher than 95% of the PLV values of the surrogate pairs, it was retained as a weight between the nodes i and j, otherwise it was set to zero. This process was executed for all four frequency bands and hence for each 20 s EEG epoch, we obtained four functional networks.

EEG Metrics
Following the computation of the undirected, weighted functional networks we characterized their structure by using brain network metrics from graph theory such as strength, eigenvector centrality and betweenness centrality [16,17]. The strength of a node is the sum of its weights whilst the betweenness centrality of a node denotes the number of times for which the shortest path lengths of other nodes pass through that node. The eigenvector centrality is a measure of the importance of a node, i.e., nodes with large eigenvector centrality are connected to other central nodes. In our study, we are interested

Functional Network
We built functional networks for each 20 sec EEG epoch by associating the network nodes to electrodes, whilst the edges between the nodes i and j were derived by the phase locking value (PLV) [19][20][21], where N is the number of samples and ∆φ i,j (t k ) is the instantaneous phase difference between the EEG signals that correspond to nodes i and j. The phases were extracted from the filtered signals using the Hilbert transform. To account for the possible effects of volume conduction [22], connections at zero time lag τ i, j = arg ∑ N k=1 e i ∆φ i, j (t k ) were discarded. Further, to account for the presence of connections that are due to chance or due to the finite length of the signals we additionally used surrogate signals. In particular, for every pair of signals, we generated 99 surrogate pairs of signals using the iterative amplitude-adjusted Fourier transform (IAAFT) with 10 iterations [23,24]. Hence, for every PLV value of the (i, j) pair of nodes we obtained a distribution of 99 PLV values from the surrogate signals. If the PLV value of the (i, j) pair was higher than 95% of the PLV values of the surrogate pairs, it was retained as a weight between the nodes i and j, otherwise it was set to zero. This process was executed for all four frequency bands and hence for each 20 s EEG epoch, we obtained four functional networks.

EEG Metrics
Following the computation of the undirected, weighted functional networks we characterized their structure by using brain network metrics from graph theory such as strength, eigenvector centrality and betweenness centrality [16,17]. The strength of a node is the sum of its weights whilst the betweenness centrality of a node denotes the number of times for which the shortest path lengths of other nodes pass through that node. The eigenvector centrality is a measure of the importance of a node, i.e., nodes with large eigenvector centrality are connected to other central nodes. In our study, we are interested in global network properties and thus we computed the mean strength S M , average betweenness centrality C B , average eigenvector centrality C E which are the average of strength, betweenness and eigenvector centrality across all nodes. In addition, we computed the variance of strength S V which is the variance of the strength across all nodes.

ECG Metrics
As is clinically standard, the acquisition of the EEG signals was performed simultaneously with an ECG acquisition and hence the selected 20-s EEG segments were accompanied by two channels of ECG data. Therefore, along with the temporal evolution of the EEG metrics we sought to investigate the temporal evolution of ECG metrics. In each 20 s epoch from the two available ECG signals, we selected the signal with the best quality. In the case where neither ECG channel showed reasonable quality, we discarded the corresponding ECG epoch from the analysis. If the amount of the available ECG epochs was less than 16 (i.e., 2/3 of the expected amount of daily epochs) we discarded that patient from this analysis. This was the case for 3 out of the 12 analyzed patients (Table 1). Similarly to the EEG acquisition, the sampling frequency of the ECG signal was either 256 Hz or 512 Hz.

HR and HRV Metrics
After downsampling the ECG signal to 256 Hz, we used the Pan-Tompkins algorithm [25] to detect the QRS complexes. Using detected R-peaks of the QRS complexes, we computed the heart rate (HR) and heart rate variability (HRV) for each ECG epoch. We estimated the HR by multiplying by three the number of R-peaks of the corresponding ECG segment (i.e., 3 × 20 s = 1 min), whilst the HRV was estimated by the root mean square of the successive difference (RMSSD) of the R-peaks [26,27].

Statistical Analysis
For each EEG and ECG metric, we obtained a distribution of values whose number was equal to the number of the analysed 20 s epochs. To account for local, single-sample fluctuations that occur in the distributions of the EEG/ECG measurements, all values were smoothed with a moving average backward filter of lag 5 (i.e., 5 samples) and normalized in the [0, 1] range. In addition, to evaluate whether the distribution of the EEG/ECG metrics of the day before the seizure occurrence was statistically significantly larger or smaller compared to the corresponding distributions of the previous days we performed multiple one-sided non-parametric Wilcoxon rank-sum tests. For instance, if a patient had six seizurefree days before the seizure occurrence, we applied five one-sided Wilcoxon rank-sum tests between the EEG metrics of the day prior to the seizure, i.e., d −1 and the EEG metrics of each previous day To decide the side of the non-parametric test (i.e., right or left sided) we compared the medians of the EEG metrics distributions from the day d −1 and the furthest day from , e.g., day d −6 . If median(d −1 ) > median(d −6 ) we applied one-sided (right tail) tests, otherwise, one-sided (left tail) tests. Results deemed significant for p-values < 0.05. To account for the effect of multiple comparisons (number of EEG/ECG metrics, number of frequency bands, number of seizure-free days prior to the seizure occurrence) we applied a Benjamini-Hochberg false discovery rate correction [28] at a significance level of 5%. In the case for which the EEG/ECG metrics of the day before the seizure were statistically significantly higher or lower compared to the corresponding metrics of all the other previous days, we performed a bootstrapping approach to ensure that the results were not due to chance. In particular, we randomly shuffled 100 times the metrics values across different days (i.e., destroying the effect of the day whilst keeping the same number of points in each day) and computed the number of times for which we observed the same pattern (i.e., the metrics of the day before the seizure to be higher or lower compared to the previous days). From this empirical distribution, we then obtained a p-value. Results were deemed significant if p − value < 0.05. All calculations and data visualization were performed in MATLAB and Statistics Toolbox Release 2020b (The MathWorks, Inc., Natick, Massachusetts, United States) and Python version 3.10.0 (Python: A dynamic, open source programming language. Python Software Foundation. URL https://www.python.org/, accessed date 9 August 2022). In addition, we used the Brain Connectivity Toolbox [17] as well as the Raincloud Plots multi-platform tool [29].

Temporal Evolution of the EEG and ECG Metrics
To gain insight into how brain network metrics vary across time and in particular whether they present any particular behavior close to seizure occurrence we started our analysis by illustrating the temporal evolution of the computed EEG metrics. Figure 2 demonstrates the temporal evolution of the brain network metrics (mean strength S M , variance of strength S V , average betweenness centrality C B , average eigenvector centrality C E ) in the alpha band, as well as the ECG metrics (HR and HRV) for the two participants (KCL1 and KCL2) who had the longest data available for analysis. We make three observations. First, there is high variability in both EEG and ECG metrics across time. This variability occurs regardless of the time of the day or the period of the seizure occurrence. Second, in the KCL 1 participant, the mean strength S M and variance of strength S V tend to take larger values in the days closer to seizure occurrence, i.e., days -1 and -2 show larger values compared to days -11 and -10. Third, all metrics and in particular the ECG metrics, show a partly periodic behavior, that can be possibly attributed to the presence of daily rhythms. ). In addition, we used the Brain Connectivity Toolbox [17] as well as the Raincloud Plots multi-platform tool [29].

Temporal Evolution of the EEG and ECG Metrics
To gain insight into how brain network metrics vary across time and in particular whether they present any particular behavior close to seizure occurrence we started our analysis by illustrating the temporal evolution of the computed EEG metrics. Figure 2 demonstrates the temporal evolution of the brain network metrics (mean strength , variance of strength , average betweenness centrality , average eigenvector centrality ) in the alpha band, as well as the ECG metrics (HR and HRV) for the two participants (KCL1 and KCL2) who had the longest data available for analysis. We make three observations. First, there is high variability in both EEG and ECG metrics across time. This variability occurs regardless of the time of the day or the period of the seizure occurrence. Second, in the KCL 1 participant, the mean strength and variance of strength tend to take larger values in the days closer to seizure occurrence, i.e., days -1 and -2 show larger values compared to days -11 and -10. Third, all metrics and in particular the ECG metrics, show a partly periodic behavior, that can be possibly attributed to the presence of daily rhythms.    delta, theta, and beta ( Figures S1-S3). This suggests that the temporal evolution of EEG and ECG metrics is highly influenced by processes that occur on various timescales with strong contributions to daily rhythms. metrics in all frequency bands. Figure 3 depicts the normalized power spectral densities of all metrics in the alpha frequency band for all study participants. We observe that the periodograms of the EEG metrics and particularly the variance of strength ( ) have in almost all participants a dominant peak around 24 h. In addition, there is periodicity in the subharmonics around 12 and 8 h. The ECG metrics show strong periodicity around 24 h across all patients. Interestingly, similar results hold for the other frequency bands, i.e., delta, theta, and beta ( Figures S1-S3). This suggests that the temporal evolution of EEG and ECG metrics is highly influenced by processes that occur on various timescales with strong contributions to daily rhythms. : variance of strength. Note that participants KCL 3, KCL 5 and KCL 12 do not have ECG metrics due to the absence of available ECG recordings (see Table 1).

Daily Distributions of the EEG and ECG Metrics
Having observed strong contributions of the daily rhythms in the temporal evolution of the EEG and ECG metrics ( Figure 3) we sought to explore the temporal evolution of the daily distributions of the EEG and ECG metrics. In particular, we investigated whether the distribution of the EEG/ECG metrics on the day prior to seizure occurrence is larger or smaller compared to the daily distributions of the corresponding metrics of all previous days. Hence, if a participant had five seizure-free days prior to seizure occurrence, i.e., { , , , , } we obtained for each metric five daily distributions and applied four one-sided Wilcoxon rank sum tests (corrected for multiple comparisons, see Section 2.6) between the metrics of the day before the seizure and the metrics of each previous day, i.e., { , , , }. We performed this analysis for all EEG/ECG metrics in all frequency bands. Figure 4 illustrates the daily distributions of the variance of strength in the alpha frequency band for all patients. We found that in 7 out of 12 participants the distribution of the day before the seizure was either statistically significant larger (KCL 1, KCL 4, KCL 10, KCL 11) or smaller (KCL 5, KCL 8, KCL 12) compared to the distributions of all previous days. Interestingly, in 2 out of 11 patients (KCL 2, KCL 9) the distributions of the two last days before the seizure occurrence were similar but statistically signifi- Figure 3. Periodograms of all the EEG and ECG metrics in the alpha frequency band. HRV: hear rate variability, HR: hear rate, C E : eigenvector centrality, C B : betweenness centrality, S M : mean strength, S V : variance of strength. Note that participants KCL 3, KCL 5 and KCL 12 do not have ECG metrics due to the absence of available ECG recordings (see Table 1).

Daily Distributions of the EEG and ECG Metrics
Having observed strong contributions of the daily rhythms in the temporal evolution of the EEG and ECG metrics ( Figure 3) we sought to explore the temporal evolution of the daily distributions of the EEG and ECG metrics. In particular, we investigated whether the distribution of the EEG/ECG metrics on the day prior to seizure occurrence is larger or smaller compared to the daily distributions of the corresponding metrics of all previous days. Hence, if a participant had five seizure-free days prior to seizure occurrence, i.e., {d −1 , d −2 , d −3 , d −4 , d −5 } we obtained for each metric five daily distributions and applied four one-sided Wilcoxon rank sum tests (corrected for multiple comparisons, see Section 2.6) between the metrics of the day before the seizure d −1 and the metrics of each previous day, i.e., {d −2 , d −3 , d −4 , d −5 }. We performed this analysis for all EEG/ECG metrics in all frequency bands. Figure 4 illustrates the daily distributions of the variance of strength S V in the alpha frequency band for all patients. We found that in 7 out of 12 participants the S V distribution of the day before the seizure was either statistically significant larger (KCL 1, KCL 4, KCL 10, KCL 11) or smaller (KCL 5, KCL 8, KCL 12) compared to the S V distributions of all previous days. Interestingly, in 2 out of 11 patients (KCL 2, KCL 9) the S V distributions of the two last days before the seizure occurrence were similar but statistically significantly smaller to all previous days. In addition, in 3 out of 12 participants the daily distributions of the S V values in the day before the seizure occurrence did not differ from the previous days. All p-values were corrected for multiple comparisons and reported in Table S1. cantly smaller to all previous days. In addition, in 3 out of 12 participants the daily distributions of the values in the day before the seizure occurrence did not differ from the previous days. All p-values were corrected for multiple comparisons and reported in Table S1. We performed this analysis for the remaining frequency bands, i.e., delta ( Figure S4), theta ( Figure S5) and beta ( Figure S6), and found that the S V distribution of the day prior to seizure occurrence was not consistently larger or smaller compared to the distributions of the previous days. Similar findings were also found for the mean strength ( Figures S7-S10), average betweenness centrality (Figures S11-S14) and eigenvector centrality (Figures S15-S18) across all frequency bands. Therefore, the only informative EEG metric that showed significant changes in the EEG metrics daily distribution prior to seizure occurrence was the variance of strength in the alpha frequency band. To also ensure that these findings were not due to chance, we performed a bootstrap method by randomly shuffling 100 times the S V values across different days (i.e., destroying the effect of the day whilst keeping the same number of points in each day). From this approach, we excluded participants whose distributions were similar across days, i.e., KCL 3, KCL6 and KCL7. All p-values < 0.05 and hence results were not attributed to chance.
Subsequently, we performed the same analysis for the ECG metrics. Figure 5 illustrates the daily distributions of the heart rate HR across the patients for whom there were available ECG recordings (Table 1). In 4 out of 9 participants, the distribution of the day before the seizure was statistically significantly larger (KCL 2, KCL7, KCL 9) or smaller (KCL 8) compared to the HR distributions of all previous days. In one out of 9 participants (KCL 4) the HR distributions of the two last days before the seizure occurrence were similar, but statistically significantly larger from the HR distribution of their previous day. For the remaining four participants either the daily HR distributions were similar across days (KCL 6, KCL 10, KCL 11) or there was a fluctuation in the HR values across days (KCL 1). The corresponding heart rate variability (HRV) distributions were qualitatively similar however with the opposite direction from the HR distributions ( Figure S19).

Lateralization of the Seizure Focus
Having investigated the temporal evolution of the daily distributions of the EEG metrics prior to seizures and found that the alpha frequency band is the most informative, we sought to explore whether we can use the EEG metrics in the alpha band for further analysis. We investigated whether the analyzed data could be informative to lateralize the hemisphere of the seizure focus. Five study participants had a clear seizure focus (Table 1; KCL 1, KCL 8 had a seizure focus on the right hemisphere; KCL 4, KCL 5, KCL 6 had the seizure focus on the left hemisphere). Figure 6 depicts the distributions of the average eigenvector centrality of each hemisphere across all analyzed patients. We found that in four out of five patients the hemisphere that contained the seizure focus had statistically significantly higher eigenvector centrality compared to the other hemisphere. We performed this analysis for the other EEG metrics, i.e., mean strength, variance of strength, betweenness centrality, however, in mean strength and variance of strength there were at least 3 out of 5 patients for which there was not a statistically significant difference between the distributions of the two hemispheres ( Figures S20 and S21). In addition, in the average between centrality in 2 out of 5 participants there were statistically significant differences between the two hemispheres and in one participant the difference was in another direction ( Figure S22).
Biomedicines 2022, 10, x FOR PEER REVIEW 10 of 16 Figure 5. Daily distributions of the heart rate HR for each seizure free day prior to seizure occurrence. Dots illustrate the HR values obtained from each analyzed ECG segment, whilst histograms and boxplots depict their distribution. Horizontal lines in the boxplots indicate the median. The day before the seizure occurrence, i.e., is denoted with blue. Days whose distributions are statistically significantly different (one-sided Wilcoxon rank sum test) from the distribution of day are illustrated in red, otherwise are depicted in green.

Figure 5.
Daily distributions of the heart rate HR for each seizure free day prior to seizure occurrence. Dots illustrate the HR values obtained from each analyzed ECG segment, whilst histograms and boxplots depict their distribution. Horizontal lines in the boxplots indicate the median. The day before the seizure occurrence, i.e., d −1 is denoted with blue. Days whose distributions are statistically significantly different (one-sided Wilcoxon rank sum test) from the distribution of day d −1 are illustrated in red, otherwise are depicted in green.
formed this analysis for the other EEG metrics, i.e., mean strength, variance of strength, betweenness centrality, however, in mean strength and variance of strength there were at least 3 out of 5 patients for which there was not a statistically significant difference between the distributions of the two hemispheres ( Figures S20 and S21). In addition, in the average between centrality in 2 out of 5 participants there were statistically significant differences between the two hemispheres and in one participant the difference was in another direction ( Figure S22).

Discussion
In this study, we examined the multiday temporal evolution of brain network metrics prior to seizure occurrence. Using samples of scalp EEG recordings from every single hour Figure 6. Distributions of the Average Eigenvector Centrality C E of each hemisphere in the alpha frequency band. Each dot depicts the C E value of the right (green) or left (orange) hemisphere as computed from a single EEG segment. KCL 1 and KCL 8 participants had the seizure focus on the right hemisphere (p-values 3.13 × 10 −6 and 3.54×, respectively, one-sided Wilcoxon rank-sum test), whilst KCL 5, KCL 6 and KCL 4 had the seizure focus on the left hemisphere (p-values 1.75 × 10 −14 , 3.63 × 10 −7 and 0.9, respectively, one-sided Wilcoxon rank-sum test).

Discussion
In this study, we examined the multiday temporal evolution of brain network metrics prior to seizure occurrence. Using samples of scalp EEG recordings from every single hour of each day we inferred functional brain networks in four frequency bands across multiple days. We quantified the topological properties of the brain networks using metrics from graph theory, i.e., mean strength, variance of strength, average betweenness centrality and average eigenvector centrality. We found that the brain network metrics fluctuate over time ( Figure 2) and exhibit periodic-like behaviour with prominent circadian features ( Figure 3 and Figures S1-S3). Considering the contribution of daily rhythms in the temporal evolution of the brain network metrics, we investigated whether the metrics in the day prior to seizure were statistically significantly larger or smaller compared to all previous days. We found that in 7 out of 12 patients the variance of strength in the alpha band on the day prior to seizure was significantly different to all previous days, and in 2 out of 12 patients the two last days prior to the seizure had similar distributions but were significantly different to all previous days (Figure 4). Having observed that the alpha frequency band is the most informative band in the temporal evolution of brain network metrics, we sought to explore whether the graph metrics can be used to lateralize the seizure focus in the alpha band. We found that in 4 out 5 patients who had a clear seizure focus the eigenvector centrality was able to lateralize the hemisphere of the seizure focus ( Figure 6).
Tools from network neuroscience have been widely used to characterize topological properties of functional networks inferred from EEG recordings [10]. In the context of epilepsy, Khambhati et al. [30] showed that long-term frequency-dependent reorganization of interictal functional networks reflects seizure proneness. In addition, Chowdhury et al. 2014 [31] showed that mean degree (equivalent to strength) in the low alpha band is able to distinguish PWE from controls whilst Pegg et al. [32] showed that mean degree can distinguish people with well-controlled and drug-resistant focal epilepsy. Furthermore, the variance of strength has been reported to be higher in PWE and their first-degree relatives compared to controls [31]. Eigenvector centrality is closely associated with computational metrics that quantify the tendency of a brain network to generate seizures and it is able to identify brain regions responsible for seizure generation [33,34]. Moreover, betweenness centrality has been reported as a useful metric to identify brain regions that are neighboring to the seizure onset zone [35]. Hence, we used in this study four network metrics, i.e., mean strength (quantifies network connectedness), variance of strength (variation in the node strength distribution), average eigenvector centrality (quantifies presence of network hubs) and average betweenness centrality (quantifies presence of nodes that act like network bridges), to capture different topological properties of the functional networks.
We found that the brain network metrics fluctuate over time and show periodic behaviour with a peak period around 24 h (Figure 2, Figure 3 and Figures S1-S3). These findings are in line with previous studies that investigated the temporal evolution of brain network metrics using intracranial EEG recordings. Kuhnert et al. [13] and Geier et al. [15] inferred functional brain networks and used strength, betweeness centrality, clustering coefficient, and shortest path length to characterize the topological properties of the networks. They reported temporal fluctuations of the metrics over time to a greater or lesser extent as well as the strong contribution of daily rhythms in this temporal variation. However, those studies did not perform their analysis on specific frequency bands but applied broadband filtering.
When we examined the temporal evolution of the daily distributions of brain network metrics, we asked whether the metric's distribution in the day prior to seizure takes smaller or larger values compared to all previous days (Figure 3 and Figures S4-S19), and found that the most informative metric was the variance of strength in the alpha frequency band (Figure 3). In 7 out of 12 participants the distribution of the variance of strength in the day prior to seizure was consistently larger (4 participants) or smaller (3 participants) compared to all previous days. In 2 out of 12 patients, the variance of strength distributions of the two last days prior to seizure was similar but statistically significantly smaller to all previous days. The variance of strength quantifies the variability of the strength values across the network nodes. High variance of strength means that the network has nodes whose strength is much higher than the average strength of the network. This might indicate a more regular network topology with nodes that act like hubs. On the other hand, low variance of strength means that the strength values across nodes are similar. Lower strength variability could be linked to increased synchronization (i.e., all nodes tend to connect together tightly and have similar large strength values) or decrease in synchronization (i.e., the connections between the nodes become looser and all nodes have similar small strength values). It is suspected that seizures are generated due to the imbalance between inhibition and excitation [36]. Previous studies that used intracranial EEG recordings and studied changes in phase synchronization and network topology metrics before, during and after seizure in a high temporal resolution (from seconds and minutes to hours) reported both a decrease and increase in synchronization prior to seizure occurrence, however, a synchronization decrease is the most prevalent [36][37][38]. In addition, Schindler at al. [12] reported that the network topology changes from a more random structure to more regular in seizures and returns back to randomness towards seizure termination.
In patients with a clear seizure focus we also found that average eigenvector centrality was higher in the hemisphere of the seizure focus ( Figure 6). Eigenvector centrality is a measure of node importance and in the context of epilepsy surgery, it has been shown to serve a useful marker for the identification of network hubs and brain regions responsible for seizure generation [33,34]. In addition, Coito et al. 2015 [39] studied directed connec-tivity in patients with left and right temporal lobe epilepsy and found different patterns of time-varying connectivity between the two groups as well as the concordance of the highest outflow region with the epileptogenic zone.
Considering that the EEG acquisition occurred simultaneously with an ECG acquisition we also investigated the multiday temporal evolution of heart rate and heart rate variability metrics. We quantified the heart rate variability using the root mean square of the successive difference of R-peaks (RMSSD). Previous studies have shown that RMSSD correlates with ECG frequency metrics, and it is a reliable metric even in ultra short-term ECG recordings [26,27]. We found that the HR and HRV have a strong periodicity of around 24 h (Figure 3). This is in line with recent work from Karoly et al. [40] and Gregg et al. [41] that used wearable devices to investigate heart rate cycles in people with epilepsy and found the presence of circadian and multiday cycles. In addition, we found that in 4 out of 9 participants the HR distribution of the day before the seizure was statistically significantly larger (3 participants) or smaller (1 participant) compared to the distributions of all previous days. In one out of 9 participants (KCL 4) the HR distributions of the two last days before the seizure occurrence was similar, but statistically significantly larger than the HR distribution of the previous day. The corresponding heart rate variability (HRV) distributions were qualitatively similar to the HR distributions however with opposite trend ( Figure S19).
Our findings should be interpreted by taking into account the study s limitations. First, the study participants were under combination therapy with two or more antiepileptic drugs (AEDs) and during their stay in the EMU there was an ongoing change in the AEDs in a patient specific manner (see Table S2). In addition, the analysed cohort was small and hence further studies are needed to extend and validate our findings.
In conclusion, the findings of this study suggest that multiday brain network metrics computed from EEG recordings could potentially be used to characterize brain network changes that occur prior to seizures. In addition, they show promise as potential biomarkers for the estimation of seizure risk and their incorporation in the study pipelines for long-term wearable EEG mobile monitoring for epilepsy forecasting and management.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biomedicines10102662/s1, The Supplementary material contains Table S1: Table with the original p-values (one-sided Wilcoxon rank-sum test) between the variance of strength distribution of the day prior to seizure (d −1 ) and the corresponding distribution of all previous days; Table S2: Table with the Patients Antiepileptic drug treatment before their admission to the hospital and the drug changes during patients stay to the hospital; Figure S1: Periodograms of all the EEG and ECG metrics in the delta frequency band; Figure S2: Same as Figure S1, but in the theta band; Figure S3: Same as Figure S1, but in the beta band; Figure S4: Daily distributions of the variance of strength S V (delta frequency band) for each seizure free day prior to seizure occurrence; Figure S5: Same as Figure S4 for the theta band; Figure S6: Same as Figure S4 for the beta band; Figure  S7: Same as Figure S4 for the mean strength and delta band; Figure S8: Same as Figure S4 for the mean strength and theta band; Figure S9: Same as Figure S4 for the mean strength and alpha band; Figure S10: Same as Figure S4 for the mean strength and beta band; Figure S11: Same as Figure S4 for the average betweenness centrality and delta band; Figure S12: Same as Figure S4 for the average betweenness centrality and theta band; Figure S13: Same as Figure S4 for the average betweenness centrality and alpha band; Figure S14: Same as Figure S4 for the average betweenness centrality and beta band; Figure S15: Same as Figure S4 for the average eigenvector centrality and delta band; Figure  S16: Same as Figure S4 for the average eigenvector centrality and theta band; Figure S17: Same as Figure S4 for the average eigenvector centrality and alpha band; Figure S18: Same as Figure S4 for the average eigenvector centrality and beta band; Figure S19: Daily distributions of the heart rate variability HRV for each seizure free day prior to seizure occurrence; Figure S20: Distributions of the Mean Strength S M of each hemisphere in the alpha frequency band; Figure S21: Distributions of the Variance of Strength S V of each hemisphere in the alpha frequency band; Figure S22: Distributions of the Betweenness Centrality C B of each hemisphere in the alpha frequency band.