Rapidly spreading seizures arise from large-scale functional brain networks in focal epilepsy

It remains unclear whether epileptogenic networks in focal epilepsy develop on physiological networks. This work aimed to explore the association between the rapid spread of ictal fast activity (IFA), a proposed biomarker for epileptogenic networks, and the functional connectivity or networks of healthy subjects. We reviewed 45 patients with focal epilepsy who underwent electrocorticographic (ECoG) recordings to identify the patients showing the rapid spread of IFA. IFA power was quantified as normalized beta-gamma band power. Using published resting-state functional magnetic resonance imaging databases, we estimated resting-state functional connectivity of healthy subjects (RSFC-HS) and resting-state networks of healthy subjects (RSNs-HS) at the locations corresponding to the patients' electrodes. We predicted the IFA power of each electrode based on RSFC-HS between electrode locations (RSFC-HS-based prediction) using a recently developed method, termed activity flow mapping. RSNs-HS were identified using seed-based and atlas-based methods. We compared IFA power with RSFC-HS-based prediction or RSNs-HS using non-parametric correlation coefficients. RSFC and seed-based RSNs of each patient (RSFC-PT and seed-based RSNs-PT) were also estimated using interictal ECoG data and compared with IFA power in the same way as RSFC-HS and seed-based RSNs-HS. Spatial autocorrelation-preserving randomization tests were performed for significance testing. Nine patients met the inclusion criteria. None of the patients had reflex seizures. Six patients showed pathological evidence of a structural etiology. In total, we analyzed 49 seizures (2-13 seizures per patient). We observed significant correlations between IFA power and RSFC-HS-based prediction, seed-based RSNs-HS, or atlas-based RSNs-HS in 28 (57.1%), 21 (42.9%), and 28 (57.1%) seizures, respectively. Thirty-two (65.3%) seizures showed a significant correlation with either seed-based or atlas-based RSNs-HS, but this ratio varied across patients: 27 (93.1%) of 29 seizures in six patients correlated with either of them. Among atlas-based RSNs-HS, correlated RSNs-HS with IFA power included the default mode, control, dorsal attention, somatomotor, and temporal-parietal networks. We could not obtain RSFC-PT and RSNs-PT in one patient due to frequent interictal epileptiform discharges. In the remaining eight patients, most of the seizures showed significant correlations between IFA power and RSFC-PT-based prediction or seed-based RSNs-PT. Our study provides evidence that the rapid spread of IFA in focal epilepsy can arise from physiological RSNs. This finding suggests an overlap between epileptogenic and functional networks, which may explain why functional networks in patients with focal epilepsy frequently disrupt.


Introduction
An increasing number of studies have provided evidence for the concept of focal epilepsy as a network disease Fig. 1. Representative rapid spread of ictal fast activity (IFA) observed in Patients 1 (upper row) and 9 (lower row). Black and yellow circles in the left column represent electrode locations on the patients' brains. The right column shows ictal electrocorticographic data of all electrodes. The traces highlighted in yellow indicate data from the top 20 electrodes ranked by IFA power, colored yellow in the left column. Data were high-pass filtered at 0.5 Hz. The x-axis represents the time relative to the seizure onset. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) , contribute to seizure genesis. However, it is unknown whether such networks develop on the brain networks subserving normal functions ( Smith and Schevon, 2016 ;Spencer et al., 2018 ). Although recent calcium imaging studies on animal models have demonstrated that epileptic and functional activities share propagation networks ( Liou et al., 2018 ;Rossi et al., 2017 ), evidence on human subjects is scarce. In epileptic patients, the association between seizures and functional networks has been evident only in a special seizure type, reflex seizure ( Koepp et al., 2016 ;Salek-Haddadi et al., 2009 ). There is an ongoing debate on whether this link could be generalized to other common types of seizures ( Avanzini et al., 2012 ;Vivekananda et al., 2019 ). Recent neuroimaging studies have found that functional network disturbance is responsible for dyscognitive and neuropsychiatric symptoms in patients with focal epilepsy ( Englot et al., 2016 ;Spencer et al., 2018 ). It has been proposed that repeated seizure discharges may alter functional networks ( Englot et al., 2016 ;Spencer et al., 2018 ); however, there is no direct evidence. Therefore, exploring the relationship between seizure propagation pathways and functional networks is vital for understanding the pathomechanisms of cognitive and neuropsychiatric dysfunction in patients with focal epilepsy.
In the last two decades, functional magnetic resonance imaging (fMRI) studies have discovered and extensively analyzed spontaneous synchronization between functionally related regions during the resting state, termed resting-state functional connectivity (RSFC) ( Power et al., 2014 ). These studies have demonstrated a close link between RSFC and brain activation during cognitive, perceptual, and motor tasks. For instance, Cole et al. reported that activity flow estimated via RSFC could predict cognitive task-induced activation ( Cole et al., 2016 ). Spatial pattern analyses of RSFC have consistently identified large-scale networks, called resting-state networks (RSNs), which correspond well to task-induced coactivation patterns ( Smith et al., 2009 ) and unique functional systems, such as visual, somatomotor, attention, and default mode ( Power et al., 2011 ;Yeo et al., 2011 ). Therefore, we hypothesized that if IFAs propagate rapidly like normal brain functional activities, they would be related to RSFC and RSNs. Although recent studies revealed that patients with focal epilepsy show strong RSFC between seizure onset and spread zones ( Bandt et al., 2014 ;Korzeniewska et al., 2014 ;Lagarde et al., 2018 ), it can reflect abnormal non-functional networks formed due to epileptic activity or underlying diseases. Hence, we compared the rapidly spreading IFAs with RSFC and RSNs of healthy subjects to test the above hypothesis.

Patients
We retrospectively reviewed data of 45 consecutive patients with medically intractable focal epilepsy who underwent electrocorticographic (ECoG) recordings as presurgical evaluation between January 2007 and December 2017 at the Kyushu University Hospital. We selected patients who showed the rapid spread of IFA ( Fig. 1 ) using the following criteria: (1) aged over 15 years, and (2) clearly visible rhythmic activity > 13 Hz spreading further than 50 mm and to more than five electrodes within 1 s ( Singh et al., 2015 ) in multiple seizures. The exclusion criteria were as follows: (1) large brain deficits or lesions, (2) apparent multiple epileptogenic foci (e.g., tubercular sclerosis), and (3) insufficient recording quality for quantitative signal analysis. For each patient, all available seizures without severe artifacts within the first 10 s were analyzed. Two board-certified epileptologists selected the patients and seizures according to the above criteria, and disagreement was resolved by discussion. This study was approved by the Kyushu University Institutional Review Board for Clinical Research.

Acquisition of ECoG
The location of electrode placement was determined based on clinical information, including seizure semiology, neurological examination, electroencephalography (EEG), magnetoencephalography, neuroimaging (such as MRI-visible lesion), and neuropsychological tests. The recording duration was based on clinical needs. ECoG data were acquired using the Nihon Kohden EEG system (Nihon Kohden, Tokyo, Japan). The sampling rate was 1000 or 2000 Hz. The reference was a scalp electrode or the average of two scalp electrodes.

Quantitative analysis of IFA
For each seizure, a board-certified epileptologist annotated the putative seizure onset time, defined as the onset of IFA at any one electrode. We excluded electrodes with artifacts from this analysis. We set the putative time of seizure onset to 0, extracted ECoG data from − 80 s to 20 s, and set the baseline period as − 80 s to − 20 s. Raw data were notch filtered. To analyze IFA in an objective and reproducible way, we computed the beta-gamma power using the following steps ( Fig. 2 ): (1) band-pass filtering at 14 Hz to 300 Hz, (2) spectral equalization via firstorder differencing to compensate for spectral roll-off ( Schindler et al., 2007 ), (3) rectification, (4) normalization by subtraction of mean value and division by standard deviation value of the baseline period, (5) estimating the time of "objective" seizure onset (see below), (6) averaging for the 4 s time duration windows centered at time points from 10 s before to 10 s after the "objective" seizure onset (1 s step), and (7) division by the maximum value across channels at the time window of 0-4 s. Before normalization, baseline data were divided into non-overlapping 1 s duration periods, and those having interictal discharges were excluded. In step (5), we applied a median filter of 200 ms to beta-gamma power, detected continuous periods longer than 3 s in which all data exceeded 0 in each channel, and defined the start of the earliest one as the "objective " seizure onset. As IFA is known to last for several seconds ( Bartolomei et al., 2008 ), we used the beta-gamma power at the time window of 0-4 s for the primary analysis. For simplification, we termed this measure as the IFA power. The other beta-gamma power data were only used to characterize the IFA time course.

Characterization of IFA time course
To describe the collective time courses of IFA, we averaged the betagamma power of the electrodes having the 10 most significant IFA powers. To illustrate the time courses of IFA spread, we applied the empirical thresholds to the beta-gamma power and treated the above-threshold values as IFA invasion. The empirical thresholds were determined as follows. An epileptologist reviewed the ictal ECoG data of all patients (the two earliest seizures in each patient) and judged whether each channel showed IFA in > 50% in the time window of 0-4 s. Then, we used a receiver operating characteristic (ROC) curve to determine the optimal cut-off value of the beta-gamma power. The values that were of 1, 1.5, and 2 times of the optimal cut-off were used as the empirical thresholds.

Electrode localization
For each patient, we localized ECoG electrodes using the method implemented in the Fieldtrip toolbox ( Stolk et al., 2018 ). We used the post-implantation computed tomography (CT) image to identify the coordinates of individual electrodes. The pre-implantation highresolution T1 weighted MRI and the post-implantation CT images were co-registered, and brain shift caused by implantation surgery was corrected ( Dykstra et al., 2012 ). We spatially normalized the T1 weighted MRI of the patients to the standard Montreal Neurological Institute (MNI) space using the Computational Anatomy Toolbox ( http://www.neuro.uni-jena.de/cat/ ) and obtained the parameters for transformation between individual and standard spaces.

Estimation of RSFC of healthy subjects (RSFC-HS)
We used resting-state fMRI (rs-fMRI) data from unrelated 100 healthy subjects from the Human Connectome Project (HCP) S500 release ( Smith et al., 2013 ;Van Essen et al., 2012 ) to estimate RSFC of healthy subjects (RSFC-HS) between electrode locations ( Fig. 2 ). Details of the HCP data acquisition and preprocessing can be found elsewhere Smith et al., 2013 ;Van Essen et al., 2012 ). We used the data provided after denoising with an independent component analysis-based method (ICA-FIX). We performed the following additional preprocessing steps to remove the artifacts more aggressively ( Smith et al., 2013 ): (1) regressing the nuisance signals, including gray matter mean signal, white matter mean signal, and cerebrospinal fluid mean signal; (2) band-pass filtering at 0.009 Hz to 0.08 Hz. We created spherical regions of interest (ROIs) with radii of 5 mm (i.e., half of the typical inter-electrode distance) centered at the electrode coordinates and transformed them into the MNI space. We calculated Pearson's correlation coefficient between the mean preprocessed rs-fMRI time series of ROIs and averaged across 100 subjects to estimate the RSFC-HS between the electrode locations.

IFA prediction based on RSFC-HS (RSFC-HS-based prediction)
We performed a recently proposed analysis named activity flow mapping ( Cole et al., 2016 ;Ito et al., 2020 ) to evaluate whether RSFC-HS between electrode locations can predict IFA power ( Fig. 2 ). For a given electrode, we estimated the activity flow from the other electrodes by RSFC and predicted the IFA power as follows: where P j is the predicted IFA power based on RSFC-HS (RSFC-HS-based prediction) of electrode j, A i is the actual IFA power of electrode i , and F ij is the RSFC-HS between the locations of electrode i and electrode j . The Brain Activity Flow Toolbox ( https://colelab.github.io/ActflowToolbox/ ) was used for this calculation. Then, the prediction accuracy was measured by comparing the actual and predicted IFA power, the statistical significance of which can be tested by methods that account for spatial autocorrelation ( Ito et al., 2020 ) (see Section 2.10 ). Although other approaches can be used to correlate RSFC with IFA ( Shah et al., 2019 ), we chose activity flow mapping analysis because it is based on computational simulation, is simple to use, is easily interpretable, and has appropriate statistical tests. However, it is difficult to test the correlation between RSFC itself and other variables because of the complex dependency between RSFC values. For example, due to spatial autocorrelation, RSFC between regions A and B is not independent of that between regions A and C if C is a neighbor of B.

Identification of RSNs of healthy subjects (RSNs-HS)
We identified RSNs of healthy subjects (RSNs-HS) in the electrode space using seed-based and atlas-based methods ( Fig. 2 ). In the seed-based method, we estimated an RSN-HS as RSFC-HS between a seed and the other electrode locations (the third column in Fig. 6 ). All electrodes were used as seeds to generate all possible seed-based RSNs-HS. In the atlas-based method, we adopted coarse and fine brain parcellations created by Yeo et al., including seven and 17 RSN-HS atlases, respectively ( Fig. 4 E) ( Yeo et al., 2011 ). Names for RSNs-HS were provided according to the previous studies in this research field ( Tang et al., 2020 ;Yeo et al., 2011 ). These parcellations were produced by rs-fMRI data collected from 1000 healthy subjects ( Yeo et al., 2011 ). We transformed the parcellation maps provided at the MNI space ( https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation_Yeo2011 ) into individual patients' spaces. We then defined each RSN-HS in the electrode space as the electrodes within 5 mm from each transformed RSN-HS atlas (yellow-colored electrodes in the fourth column in Fig. 5 ). Among all 24 atlas-based RSNs-HS, those with less than five electrodes were excluded from further analysis because we included IFAs invading at least five electrodes in this study.

Estimation of RSFC of patients (RSFC-PT), IFA prediction based on RSFC-PT (RSFC-PT-based prediction), and estimation of RSNs of patients (RSNs-PT)
Our primary purpose was to compare IFA with rs-fMRI databases of healthy subjects. However, we also estimated the RSFC and RSNs of each patient (RSFC-PT and RSNs-PT, respectively) and analyzed their association with IFA to confirm previous studies. Due to the retrospective nature of this study, rs-fMRI data of the patients were unavailable. Therefore, we calculated the RSFC-PT from ECoG data using an established method ( Foster et al., 2015a ;Keller et al., 2013 ). We selected a 300 s long interictal ECoG data while lying awake with minimal contamination of undesirable events, such as interictal epileptiform discharges (IEDs), pathological high-frequency oscillations, and electromyographic signals as possible. We determined the sleep stage by simultaneously recording the scalp EEG. Two board-certified epileptologists annotated undesirable events, and disagreements were resolved by discussion. The data were notch filtered, re-referenced to a common average reference, and filtered using multiple band-pass filters of 10 Hz width between 50 Hz and 150 Hz (e.g., 50-60, 60-70…140-150 Hz). Here, we focused on the high-gamma range because several ECoG studies have consis- tently reported that connectivity in this range corresponds well to that observed using rs-fMRI ( Foster et al., 2015b ;Keller et al., 2013 ). For each 10 Hz band, the signal envelope was calculated using the Hilbert transform and normalized to its mean. Normalized envelopes were averaged across 10-Hz width bands to obtain high-gamma power (HGP). We low-pass filtered the HGP signal at 1 Hz and discarded data segments from 3 s before to 3 s after any undesirable events. If the remaining data were shorter than 200 s, we excluded the patient from further analysis. Then, we calculated RSFC-PT as Pearson's correlation coefficient between low-pass filtered HGP signals. IFA prediction based on RSFC-PT (RSFC-PT-based prediction) and estimation of seed-based RSNs of each patient (RSNs-PT) were performed in the same way as for RSFC-HS-based prediction and seed-based RSNs-HS estimation, respectively, where RSFC-PT was used instead of RSFC-HS.

Statistical analysis
We evaluated the association between IFA power and RSFC-HSbased prediction, seed-based RSNs-HS, RSFC-PT-based prediction, or seed-based RSNs-PT using Spearman's rank correlation coefficients ( r s ) ( Ito et al., 2020 ; Vos de Wael et al., 2020 ) ( Fig. 2 ). We used the rank-biserial correlation coefficients ( r rb ) to compare IFA power with atlasbased RSNs-HS ( Fig. 2 ). In comparison with seed-based or atlas-based RSNs (i.e. seed-based RSNs-HS, seed-based RSNs-PT, and atlas-based RSNs-HS), we correlated IFA power with all RSNs. We then applied the maximum r s or r rb values as the test statistics because we were interested in whether IFA power corresponds to at least one RSN ( Brookes et al., 2011 ). Seizure discharges recorded by ECoG are known to have spatial autocorrelation ( Martinet et al., 2015 ), which causes inflated type I error in any statistical tests assuming independence ( Alexander-Bloch et al., 2018 ). Hence, we performed the spatial autocorrelation-preserving randomization test for the correlation measures using the BrainSpace toolbox ( Vos de Wael et al., 2020 ) ( Fig. 2 ). We randomized IFA power 5000 times using the Moran spectral randomization method ( Wagner and Dray, 2015 ) ( Fig. 2 ) and calculated the surrogate statistics. We then estimated the one-sided P value as the proportion of the surrogate values exceeding the observed value and considered a P value < 0.025 as significant. For the comparison of IFA power and seed-based or atlas-based RSNs, we retained only the maximum correlation value as a surrogate in each randomization to correct for multiple comparisons ( Brookes et al., 2011 ). For each patient, the individual seizure P values were combined into the harmonic mean P value with equal weights ( Wilson, 2019 ) to perform the patient-level analysis. For the group-level analysis, as the number of seizures varied across patients, we collected the P values of the N earliest seizures of each patient to calculate the harmonic mean P value with equal weights, where N denotes the smallest number of seizures per patient. We considered a harmonic mean P value < 0.025 to be significant in the patient-and group-level analyses.

Visualization
For effective visualization, we mapped the spatial distribution of IFA power, RSFC-based prediction, and seed-based RSNs on the patient's cortical mesh using the interpolation method implemented in the Fieldtrip toolbox ( Stolk et al., 2018 ).

Data and code availability
The patient data that support the findings of this study are not publicly available because of patient privacy but may be obtained from the authors upon reasonable request. The HCP rs-fMRI data ( https://db.humanconnectome.org/ ) and RSN-HS atlases ( https:// surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation_Yeo2011 ) are publicly available. All analyses were performed using Matlab 2017b (MathWorks, Natick, MA) along with freely available toolboxes.

Patient characteristics
Among the 45 patients, nine patients met the inclusion criteria (five males and four females; age at surgery, 16-33 years; age at onset, 2-31 years). Table 1 shows the patient characteristics. None of the patients had reflex seizures. The seizure onset electrode was located in the frontal lobe in six patients, the temporal lobe in two patients, and the parietal lobe in one patient. The MRI-visible lesions were observed in four patients (Patients 1, 2, 5, and 6), located in the parietal, temporal, frontal, and frontal lobes, respectively ( Fig. 3 A). Resective surgery following ECoG recording was performed in eight patients; the pathological findings were focal cortical dysplasia in six patients and gliosis in two patients. One patient with gliosis (Patient 8) was diagnosed with post-encephalitis epilepsy. The number of implanted electrodes was between 69 and 122. Fig. 3 B and C show the coverage of the implanted electrodes in the MNI space of all patients. The follow-up period after surgery was at least 41 months. After surgery, five patients were free of seizures (Engel class I).

Seizure characteristics
We analyzed 49 seizures that showed rapidly spreading IFA in total. The number of seizures per patient ranged from two to 13 ( Table 1 ). Therefore, for the group-level statistical tests, we collected two earliest seizures from each patient . Fig 4 A shows the spatiotemporal patterns of beta-gamma power, reflecting IFA, in the representative seizures for each patient. In Patients 1, 2, 4, 6, and 8, IFAs arose from multiple localized foci and involved relatively limited electrodes ( Fig. 4 A). In Patients 5 and 9, IFAs occurred almost simultaneously in broad regions and rapidly spread to a large portion of the electrodes ( Fig. 4 A). In Patients 3 and 7, a single region initially produced IFAs, and they propagated to the limited electrodes ( Fig. 4 A). As presented in Fig. 4 B, beta-gamma power collectively peaked within 10 s after seizure onset. The area under the ROC curve of beta-gamma power, used to predict manually identified IFA, was 0.94, and the optimal cut-off was 0.13, with a sensitivity of 83% and specificity 90%. We adopted 0.13, 0.20, and 0.26 (1, 1.5, and 2 times of the optimal threshold) as the empirical thresholds to visualize the seizure spreading characteristics. Fig. 4 C displays the mean number of electrodes with IFA estimated by the empirical thresholds to the beta-gamma power. Fig. 4 D shows the mean seizure diameter (the maximum distance between electrodes with IFA). On average, the seizure diameter rapidly increased just after seizure onset, illustrating Here, we treated electrodes as having IFA when their normalized beta-gamma power exceeded the empirical thresholds (0.13, 0.20, and 0.26). (D) Maximum distance between electrodes with IFA (seizure diameter). The x-axis represents time relative to the seizure onset. Data from 10 s before to 10 s after the seizure onset are shown in (B)-(D). Pt, patient; Sz, seizure; th, threshold. the spatially non-contiguous spreading patterns of IFAs in the analyzed seizures ( Fig. 4 D).

IFA power vs. RSFC-HS-based prediction and RSNs-HS
Between IFA power and RSFC-HS-based prediction, the mean (SD) r s value was 0.39 (0.16). At the individual seizure level, 28 (57.1%) seizures of seven patients showed significant r s values ( Fig. 5 A). In the patient-level analysis, six (66.7%) patients reached statistical sig-nificance ( Table 2 ). This correlation was significant at the group level ( Table 2 ). The second column in Fig. 6 displays RSFC-HS-based prediction in the representative seizures, where asterisks denote significant correlations with IFA power.
In comparing IFA power and seed-based RSNs-HS, the mean (SD) maximum r s value was 0.45 (0.14). Maximum r s was significant in 21 (42.9%) seizures at the seizure level ( Fig. 5 B), significant in six (66.7%) patients at the patient level, and significant at the group level ( Table 2 ). Between IFA power and atlas-based RSNs-HS, the mean (SD) maximum . Asterisks indicate significance in spatial autocorrelation-preserving randomization tests ( * P < 0.025). (D) Correlation coefficients between IFA power and each atlas-based RSN-HS averaged across seizures for each patient. The colors of the rectangles beside the RSN names refer to those in E. Asterisks indicate significance in spatial autocorrelation-preserving randomization tests ( * P < 0.025) in at least one seizure. (E) Brain maps of RSN-HS atlases ( Yeo et al., 2011 ). Control, control network; Default, default mode network, DorsAttn, dorsal attention network; Limbic, limbic network; NA, not analyzed; Pt, patient; SomMot, somatomotor network; Sz, seizure; TempPar, temporal-parietal network; VentAttn, ventral attention network; Visual, visual network. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) r rb value was 0.43 (0.12). Maximum r rb was significant in 28 (57.1%) seizures ( Fig. 5 C). In these seizures, the correlated RSNs-HS were the default mode network (DMN) and its subnetwork, the temporal-parietal network, the somatomotor network, the subnetwork of the dorsal attention network, and the control network ( Fig. 5 D). Maximum r rb was significant in seven (77.8%) patients at the patient level and significant at the group level ( Table 2 ). Taken together, 32 (65.3%) seizures showed a significant correlation with either seed-based or atlas-based RSNs-HS, but this ratio varied across patients. Twenty-seven (93.1%) of 29 seizures in six patients (Patients 1, 2, 3, 5, 6, and 9) significantly correlated with seed-based or atlas-based RSNs-HS ( Fig. 5 B-C). The third and fourth columns in Fig. 6 show seed-based and atlas-based RSNs-HS, respectively, in the representative seizures. Asterisks indicate significant correlations with IFA power.

IFA power vs. RSFC-PT-based prediction and RSNs-PT
In Patient 1, we could not extract IED-free interictal ECoG data longer than the predefined duration (200 s) for estimating RSFC-PT and RSNs-PT. Therefore, we compared the IFA power of Patient 1 with only RSFC-HS and RSNs-HS, but not RSFC-PT and RSNs-PT.
The mean (SD) r s value between IFA power and RSFC-PT-based prediction was 0.48 (0.12). Significant r s values were observed in 43 (93.3%) seizures in eight patients ( Table 2 ). This correlation was significant in eight (100%) patients at the patient level and significant at the group level ( Table 2 ).
Between IFA power and seed-based RSNs-PT, the mean (SD) maximum r s value was 0.55 (0.12). Maximum r s was significant in 38 (84.4%) Fig. 6. Visualization of representative correspondence between ictal fast activity (IFA) power of patients with focal epilepsy and resting-state functional connectivity of healthy subjects (RSFC-HS) or resting-state networks of healthy subjects (RSNs-HS). Circles on the brains represent electrode locations. The first column shows IFA power; the second column shows RSFC-HS-based prediction; the third and fourth columns show seed-based and atlas-based RSN-HS, respectively, that correlated the most with IFA power. Asterisks indicate significance in spatial autocorrelation-preserving randomization tests ( * P < 0.025). Each row represents an individual seizure. Pt, patient; RSFC-HS, resting-state functional connectivity of healthy subject; RSN-HS, resting-state network of healthy subjects; Sz, seizure. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2
Results of correlation analysis between ictal fast activity and resting-state functional connectivity or networks. seizures at the seizure level, significant in six (75.0%) patients at the patient level, and significant at the group level ( Table 2 ).

Discussion
The main new findings of our study are that there are rapidly spreading IFAs associated with RSFC-HS and RSNs-HS in patients with focal epilepsy. Furthermore, to the best of our knowledge, this is the first report describing seizure propagations along major RSNs, such as the DMN. As the rapid spread of IFA is a proposed biomarker for epileptogenic networks , our results suggest an overlap between epileptogenic and functional networks. Significant correlations with RSFC-HS or RSNs-HS were observed in two-thirds of seizures in total, but this ratio varied across patients. Such correlations were seen in most seizures in six patients and some seizures in one patient, but no such correlations were seen in any seizures in two patients. Although it is difficult to define the reason for this variability due to the small sample size, the two patients without significant correlations did not have evidence of structural etiology, whereas six of the remaining seven patients did. Therefore, a structural etiology may predispose the rapidly spreading seizures to healthy RSFC or RSNs.
We observed that the predicted IFA by activity flow mapping with the RSFC-HS was correlated with the observed one. Activity flow mapping is based on the brain network simulation and predicts local brain activities assuming that they spread through functional connections estimated by RSFC ( Cole et al., 2016 ;Ito et al., 2020 ). This analysis has uncovered the link between cognitive task-induced activations and RSFC ( Cole et al., 2016 ). Therefore, we can interpret our results that the rapid spread of IFA analyzed in this study follow RSFC-HS, such as functional brain activities. This interpretation is in line with recent calcium imaging studies on animal models, demonstrating that epileptic discharges and normal brain activities share propagation pathways ( Liou et al., 2018 ;Rossi et al., 2017 ). A growing number of studies on human patients revealed that the spread of focal seizures correlates with the structural or functional connectivity measured by the patients' MRI or electrophysiological data ( Bandt et al., 2014 ;Korzeniewska et al., 2014 ;Lagarde et al., 2018 ;Proix et al., 2017 ;Shah et al., 2019 ). Our analysis using the RSFC of each patient (RSFC-PT) confirmed these previous findings. Although the origin of these seizure propagation-related connections is currently unclear, our findings from the comparison with the RSFC-HS suggest that they are pre-existing and physiological, rather than disease-induced.
To date, the relationship between seizures and functional systems is only evident in reflex seizures ( Koepp et al., 2016 ). For instance, reading-induced seizures are known to activate brain areas physiologically subserving language and facial motor tasks ( Salek-Haddadi et al., 2009 ). However, Vivekananda et al. (2019) recently demonstrated that memory tasks increased the epileptic discharges in the memoryrelated networks in patients with focal epilepsy without reflex seizures ( Vivekananda et al., 2019 ). During the last decade, the new concept of "system epilepsy" has been proposed, hypothesizing that some types of epilepsies reflect the pathological expression of a specific brain system subserving normal physiologic functions ( Avanzini et al., 2012 ). Cumulative evidence has suggested that two types of genetic generalized epilepsies, namely absence epilepsy and juvenile myoclonic epilepsy, can be understood as system epilepsy ( Avanzini et al., 2012 ). In our study, none of the patients reported reflex seizures, and six of them showed evidence for structural etiology. Therefore, our results extend these previous studies, supposing that the interplay between epileptic and functional networks is a more common phenomenon than previously expected.
Many rs-fMRI studies on epileptic patients have consistently exhibited the disturbance of RSNs ( Englot et al., 2016 ), which may be responsible for cognitive impairment. For example, there is a close relationship between the disconnection of DMN and memory impairment in temporal lobe epilepsy ( McCormick et al., 2014 ). It has been proposed that seizure discharges may cause RSN disturbance ( Englot et al., 2016 ;Spencer et al., 2018 ). However, to the best of our knowledge, no studies before this one have clarified seizure propagations along RSNs. Our findings indicate that passing seizure discharges can induce neural damage to RSNs. It is also possible that due to repeated seizures, the RSNs of our patients became excitable as a whole to rapidly propagate seizures. Such RSN-wide excitability could lead to desynchronization because the balance of excitation and inhibition is thought to be crucial for interregional functional connectivity ( Deco et al., 2014 ). Although there is no direct evidence for such seizure-induced network-wide excitability, transcranial magnetic stimulation studies have found widespread cortical excitability in patients with focal epilepsy, which is associated with recurrent seizures ( Badawy et al., 2013 ;Hamer, 2005 ). From the viewpoint of intractability, this notion of RSN-wide excitability may explain why the rapid spread of IFA is a marker of poor surgical outcomes. In this study, three patients (Patients 1, 4, and 9) showed poor surgical outcomes. Their seizures were associated with the DMN or control network, the two largest RSNs in the atlases. In contrast, none of the five patients with good surgical outcomes showed correlations with these two RSNs. Hence, the excitability of large RSNs may worsen the outcomes of epilepsy surgery.
In pre-epilepsy surgery evaluation, intracranial electrodes are often implanted according to expected propagation networks ( Gonzalez-Martinez et al., 2014 ). Currently, these networks are mostly based on accumulated knowledge from past patients. For instance, Spencer (2002) proposed three large-scale epilepsy networks based on the extensive observation of patients with focal epilepsy: limbictemporofrontal network, occipitotemporal network, and parietofrontal network ( Spencer, 2002 ). Although it is challenging to compare directly, our study indicates that these networks may overlap with RSNs. Furthermore, the results of this study suggest that the RSFC-HS and RSNs-HS can guide seizure spreading pathways and make electrode implantation more effective. In particular, the standard brain RSN maps we used are readily available, and our method for using them can be quickly adopted and tested in a clinical setting. Recent studies on patients with intractable focal epilepsy reported that rs-fMRI is a useful presurgical evaluation ( Tracy and Doucet, 2015 ). Our results claim that rs-fMRI data from healthy subjects are also valuable in surgical planning.
This study has several limitations. First, a significant limitation of this study is the small sample size. Therefore, our findings need to be confirmed through more extensive studies. One reason for the limited sample size was that we focused on the highly selected seizure pattern. There are various seizure propagation patterns ( Kutsy et al., 1999 ;Schiller et al., 1998 ), indicating multiple coexistence mechanisms. Not only synaptic but also non-synaptic mechanisms are proposed ( Zhang et al., 2014 ). Therefore, we believe that separately analyzing each pattern is vital to explore the pathophysiology. Similar to our strategy, Martinet et al. (2015) focused on the slow recruitment in focal seizure generalization ( Martinet et al., 2015 ). Further studies are needed to test whether the findings on specific patterns can be generalized to others. Second, because of the limited number of subjects, the relationship between the affected RSN and seizure semiology was unclear. In this study, the IFA of three patients displayed correlations with DMN. Two of them had the same seizure type, focal impaired awareness automatism seizure; however, this phenotype is not specific and is shared in many patients with focal epilepsy. As it is of great interest whether there are specific seizure types for affected RSNs, larger studies are needed. Third, rs-fMRI data of the patients were unavailable. Instead, we estimated the RSFC of the patients by using the HGP of ECoG data. Several studies have confirmed that ECoG-and rs-fMRI-based RSFC are tightly linked ( Foster et al., 2015a ;Keller et al., 2013 ). Further, ECoG-based RSFC has an advantage over rs-fMRI-based RSFC in that it can eliminate the possible effects of IEDs ( Khoo et al., 2017 ), which are undetectable on the scalp. Thus, we believe that the RSFC estimation of our patients was not inferior to that based on rs-fMRI data. Currently, rs-fMRI databases of epileptic patients are being developed. In the future, these data may uncover the difference between the rs-fMRI data of patients and healthy subjects on the detectability of seizure spreading patterns in a clinical setting. Fourth, the coverage of ECoG is limited. Hence, we were blinded to the seizure discharges outside the electrode sheets and deep inside the brain. Therefore, it is possible that IFA involving broad regions from the beginning, observed in Patients 5 and 9, might represent the later stage of seizures arising from outside recording areas or deep structures such as the bottom of the sulcus. Furthermore, owing to limited coverage, we could only capture the propagation along with a part of the RSNs. It is unclear whether IFAs spread to involve all RSNs.

Conclusion
Our study provides evidence that IFA can rapidly propagate along physiological RSFC and RSNs in patients with focal epilepsy without reflex seizures. These findings suggest that epileptogenic and functional networks share underlying structures in non-reflex focal epilepsy. Our results may explain why functional networks of patients with focal epilepsy frequently disrupt.

Data and code availability statement
The patient data that support the findings of this study are not publicly available because of patient privacy but may be obtained from the authors upon reasonable request. The HCP rs-fMRI data ( https://db.humanconnectome.org/ ) and RSN atlases ( https://surfer. nmr.mgh.harvard.edu/fswiki/CorticalParcellation_Yeo2011 ) are publicly available. All analyses were performed using Matlab 2017b (Math-Works, Natick, MA) along with freely available toolboxes.