Recurrence Quantification Analysis of Dynamic Brain Networks

Evidence suggests that brain network dynamics is a key determinant of brain function and dysfunction. Here we propose a new framework to assess the dynamics of brain networks based on recurrence analysis. Our framework uses recurrence plots and recurrence quantification analysis to characterize dynamic networks. For resting-state magnetoencephalographic dynamic functional networks (dFNs), we have found that functional networks recur more quickly in people with epilepsy than healthy controls. This suggests that recurrence of dFNs may be used as a biomarker of epilepsy. For stereo electroencephalography data, we have found that dFNs involved in epileptic seizures emerge before seizure onset, and recurrence analysis allows us to detect seizures. We further observe distinct dFNs before and after seizures, which may inform neurostimulation strategies to prevent seizures. Our framework can also be used for understanding dFNs in brain function in health and other neurological disorders.


Introduction
The brain is a complex dynamic system. To map, model and study brain structure and function, it is useful to define brain networks (Fornito et al., 2016;Bassett and Sporns, 2017). The study of brain networks has transformed our understanding of the brain, and it has the potential to revolutionize the clinical management of neurological disorders (Stam, 2014). Two main types of brain networks have been considered: structural and functional networks (Bullmore and Sporns, 2009;Sporns 2013). Structural networks describe the anatomical connectivity of the brain and are relatively stable on short time scales (i.e., seconds to minutes). In contrast, functional networks are inferred from statistical dependencies between neural signals recorded from different brain regions. Statistical dependencies are then assumed to represent functional couplings between brain regions. The statistical dependencies between signals are not stationary, making functional networks time-dependent on short time scales (tens or hundreds of milliseconds) (Sporns, 2013). A growing body of evidence shows that dynamic functional networks (dFNs) capture crucial aspects underlying normal function and dysfunction of the brain (Hutchison et al., 2013;Calhoun et al., 2014;Cohen, 2018). In particular, epilepsy, which we will focus in the present study, has been considered to be a dynamical disease of the brain (da Silva, 2003), and dFNs have been useful for characterizing the epileptic brain (Lehnertz et al., 2014).
A number of different approaches have been employed to study dFNs (Braun et al., 2018). A common approach has been to calculate some measures from the functional networks and track their changes over time (Schindler et al., 2008;Kramer et al., 2010;Lehnertz et al., 2014;Fuertinger et al., 2016). For example, a time-dependent analysis of the average shortest path length and the clustering coefficient revealed that functional networks evolve from a more random topology before seizures towards a more regular topology during seizures and back to a more random topology after seizure offset (Schindler et al., 2008). This approach is limited by an a priori choice of measures that may or may not fully characterize the dynamics of the functional networks. Another common approach is to use a Bayesian framework to characterize dFNs. In particular, hidden Markov models have been employed to analyze dFNs (Eavani et al., 2013;Sourty et al., 2016;Vidaurre et al., 2018). For example, product hidden Markov models have been used to identify brain networks involved in dementia with Lewy bodies (Sourty et al., 2016). This approach makes the assumption that the dynamics of the state is Markovian, i.e. the transition between different states is a memoryless stochastic process. However, long-term correlations in temporal patterns of brain activity suggest that this assumption may not always hold (Kitzbichler et al., 2009;Chialvo, 2010;Ezaki et al., 2019). Thus, these approaches make assumptions that may hinder a comprehensive assessment of the dynamics of functional networks. Figure 1: Scheme of the data analysis procedure to apply recurrence quantification analysis to dynamic functional brain networks. (a) We segment brain activity into windows. (b) From each window, we infer a functional brain network. (c) We employ a distance measure for assessing the similarity between functional networks at each pair of time windows, resulting in a distance matrix. (d) We obtain a recurrence plot by assessing the distances between functional networks. Black dots correspond to entries in the distance matrix shown in (c) whose value is smaller than a threshold. (e) We then use recurrence quantification analysis to interrogate the recurrence plot and consequently characterize the dynamics of the functional networks. Note that, in our analysis, the windows represented in (a) overlapped with adjacent windows.

Data acquisition and pre-processing
We study two data sets from patients with epilepsy: one comprising MEG recordings and the other containing stereo EEG (sEEG) recordings.

MEG recordings
We consider resting-state MEG data recorded from 26 people with JME and 26 healthy controls. The people with epilepsy were recruited from a specialist clinic for epilepsy at University Hospital of Wales in Cardiff, and the healthy controls were volunteers who had no history of significant neurological or psychiatric disorders. The two groups were age and gender matched (age range [19,45], median 27 years, and 8 males in the epilepsy group; age range [18,48], median 27, and 7 males in the control group). People with epilepsy had a number of different seizure types and were taking anti-epileptic drugs (see Krzemiński et al., 2019 for more details about this dataset). MEG data were acquired using a 275-channel CTF radial gradiometer system (CTF System, Canada) at a sampling rate of 600 Hz. Recording sessions lasted for approximately 5 minutes per individual. The participants were instructed to sit steadily in the MEG chair with their eyes focused on a red dot on a grey background. The participants also underwent a whole-brain T1-weighted MRI acquired using a General Electric HDx 3T MRI scanner and an 8-channel receiver head coil (GE Healthcare, Waukesha, WI) with an axial 3D fast spoiled gradient recalled sequence (echo time 3 ms; repetition time 8 ms; inversion time 450 ms; flip angle 20º; acquisition matrix 256×192×172; voxel size 1×1×1 mm). This study was approved by the South East Wales NHS ethics committee, Cardiff and Vale Research and Development committees, and Cardiff University School of Psychology Research Ethics Committee. Written informed consent was obtained from all participants.

sEEG recordings
We also used a data set comprising 10 people with drug-resistant focal epilepsy who underwent invasive monitoring with stereo EEG at the 999 Brain Hospital, China. Stereo EEG is an advanced procedure in the epilepsy surgery evaluation, to help delineate the irritative and seizure onset zones, and hence decide the suitability and plan epilepsy surgery (Duncan et al., 2016). The age range of the group was [16,31], median 23, and 9 individuals were males. Electrode implantation locations were personalized according to imaging and non-invasive EEG data; the number of electrodes per implantation ranged from 5 to 16, and each electrode had 2 to 16 contacts (i.e., channels). Four individuals received bilateral implantations, two individuals had electrodes implanted in their right hemispheres, and the other four individuals received implantations in their left hemispheres. Stereo EEG data was acquired using Nihon Kohden recording system at a sampling rate of 1 or 2 kHz. For all individuals, high-resolution T1-weighted MRI were acquired before electrode implantation, and computed tomography (CT) were acquired after electrode implantation. Co-registration of the CT to the MR images allowed us to determine whether contacts were placed in grey matter. All individuals had at least 2 seizures recorded. Seizure onset was defined by epileptologists at the 999 Brain Hospital and corroborated by one of the authors (KH) who also marked seizure offset. This study was approved by 999 Brain Hospital ethics committee and South China Normal University ethics committee. Written informed consent was obtained from all participants.
We restricted our analysis to artefact-free sEEG channels placed on grey matter as established from co-registration of the CT scans to the MRIs. (Artefact-free channels were identified by KH.) The number of channels ch considered per individual ranged from 19 to 83 (median 63.5). We selected peri-ictal epochs of data containing 300 s before seizure onset (pre-ictal), seizure (ictal), and 300 s after seizure offset (post-ictal). We then neglected peri-ictal epochs whose pre-ictal 300 s overlapped with the post-ictal 300 s of the previous peri-ictal epoch. Thus, we end up with 2 to 4 peri-ictal epochs per individual (3 individuals had 2 epochs each; 5 individuals had 3 epochs each; and 2 individuals had 4 epochs each, making a total of 29 peri-ictal epochs). The data was re-referenced to the average of all artefact-free channels, filtered in a broad frequency band (0.5-120 Hz), which encapsulates the traditional clinical frequency bands (delta, theta, alpha, beta, and gamma (Buzsaki, 2006)), notch filtered to remove power line interference (48 to 52 Hz), and down-sampled to 250 Hz.

Inferring functional networks
We then constructed functional networks from the pre-processed MEG and sEEG data.

Dynamic MEG functional networks in the source space
To compute functional networks from the MEG data, we first transformed the MEG recordings from the sensor space to the source space. This procedure consisted of co-registering the MEG sensors with the structural MRI using the locations of the fiducial coils in the CTF software (MRIViewer and MRIConverter). Then, we inferred a volume conduction model from the MRI scan using a semirealistic model (Nolte, 2003). Finally, we reconstructed the source signals from the sensor signals using a linear constrained minimum variance (LCMV) beamformer on a 6-mm template with a localspheres forward model in Fieldtrip (Oostenveld et al., 2011;http://www.ru.nl/neuroimaging/fieldtrip). Sources were mapped into the 90 brain regions of the Automated Anatomical Label (AAL) atlas (Hipp et al., 2012). More details about these methods were provided in our previous study (Krzemiński et al., 2019).
To compute dFNs, we divided the 90-dimensional source reconstructed MEG signals into segments, each of which was composed of 500 samples (i.e., 2 s). Each segment was subsequently used for constructing one functional network. The choice of segment length balances the need of a sufficient number of samples to infer a reliable functional network and the need of a sufficiently large number of functional networks, , for analyzing their recurrence over time. The segment size of 500 is within a typical range in both MEG and EEG studies of functional connectivity (see e.g. Khambhati et al., 2015;Colclough et al., 2016;Stahn and Lehnertz, 2017). We set the overlap between consecutive segments to be 80% such that consecutive segments shared 400 samples. This choice represents a compromise between the need of sufficiently many networks for subsequent analysis, which is satisfied with a large overlap, and the need of avoiding trivial recurrences between consecutive functional networks. See Supplementary Material S1 for computational results underlying the choice of the 80% overlap. The MEG data had different lengths for different participants. Therefore, we only considered the first = 506 segments, which was the minimum number of segments among all participants, to exclude the potential impact of different recording lengths on our results.
For each segment, we built two functional networks using two established methods (Colclough et al., 2016): phase lag index (PLI) (Stam et al., 2007) and amplitude envelope correlation (AEC) with orthogonalized signals (Hipp et al., 2012) (see Supplementary Material S2 for more details). We selected these two methods to capture potentially complementary signatures of dFNs. Note that since we considered 4 frequency bands for each definition of functional connectivity (i.e., PLI and AEC), we obtained eight sequences of matrices ( ) = ( ( )), where , = 1, … , 90, and = 1, … , 506. Each matrix ( ) represented a functional network for segment , and matrix entry ( ) represented the strength of the functional connectivity between regions and . Each matrix ( ) was symmetric, i.e., the functional networks were undirected.

Dynamic sEEG functional networks in the sensor space
Following the same procedure as the one employed for the MEG data, we divided each of the 29 sEEG peri-ictal epochs into segments of 500 samples (i.e., 2 s) using an overlap of 400 samples between consecutive segments. Because different peri-ictal epochs contained seizures of different lengths (from 9 to 181 s, median 75 s), different peri-ictal epochs resulted in different numbers of segments (from 1542 to 1822, median 1632).
Because the sEEG data was recorded intracranially close to the brain sources, we computed functional networks in the sensor space, i.e. using the channels as network nodes and functional connections as the statistical dependencies between the signals recorded at the channels. For each segment, we inferred a functional network using the absolute value of the Pearson's correlation coefficient between pairs of channels (Rummel et al., 2010;Rummel et al., 2015;Lopes et al., 2017). Due to volume conduction, the Pearson's correlation is not a reliable measure of functional connectivity for non-invasive data such as scalp EEG or MEG data (Bastos and Schoffelen, 2016). In contrast, volume conduction is not an issue in invasive data such as sEEG, and the Pearson's correlation is a reliable method (Rummel et al., 2010), hence our choice. Thus, we obtained 29 timevarying matrices ( ) = ( ( )) of size ( ch × ch ), where = 1, … , , representing functional networks through pre-ictal, ictal, and post-ictal periods. Note that the number of channels, ch , was fixed for each individual, but varied from one peri-ictal epoch to another even for a single individual due to seizure of different lengths.

Recurrence plots and distance measures
We used recurrence plots (RPs) to characterize dFNs. For a dynamical system characterized by a vector time series ⃗( ), where = 1, … , , an × recurrence matrix is defined as where ( ⃗( 1 ), ⃗( 2 )) is a distance measure between ⃗( 1 ) and ⃗( 2 ), and is a small distance which defines an upper limit of discrepancy between recurrent states (Marwan et al., 2007;Marwan et al., 2009). The recurrence matrix is a symmetric matrix, and its main diagonal entries are equal to 1.
To compute RPs of dFNs, we replaced the vectors ⃗( ) by matrices ( ) and used distance measures for pairs of weighted networks (i.e., matrices). The recurrence matrix for a dFN is given by where ( ( 1 ), ( 2 )) is the distance between functional networks ( 1 ) and ( 2 ) measured according to a distance measure for networks. There is a variety of distance measures for networks (Livi and Rizzi, 2013), but a good choice for functional networks is unknown. We therefore used the following six distance measures.
First, we used the Frobenius norm of the difference between a pair of networks, which is given by (Kurmukov et al., 2016) (
(3) for evaluating symmetric positive-definite matrices (Arsigny et al., 2006), which our functional networks are. Third, we used the spectral norm of the difference between a pair of networks, which is given by where max {•} is the largest eigenvalue of the matrix, and * ( ) is the conjugate transpose of ( ) (Miller et al., 2012).
The other three distance measures were based on spectral properties of the Laplacian networks, specifically their Fiedler vectors. The Fiedler vector of a network is the eigenvector corresponding to the smallest positive eigenvalue of the Laplacian matrix of the network, which is often referred to as the algebraic connectivity of the network. The use of the Laplacian matrix is motivated by evidence showing that the Laplacian matrix is better for computing spectral distances than the adjacency matrix (Wilson and Zhu, 2008). The Fiedler vector characterizes the partitioning of the network into communities (Newman, 2006). Here, we considered the so-called symmetric normalized Laplacian matrix given by ′( ) = −1/2 ( ) ( ) −1/2 ( ), where ( ) = ( ) − ( ) is the combinatorial Laplacian, and ( ) is a diagonal matrix whose main diagonal entries are given by ( ) = ∑ ( ) =1 . We denote the normalized Fiedler vector of ′( ) by ⃗( ) = ( 1 ( ), 2 ( ), … , ( )) ⊤ , where ⊤ represents the transposition. We considered the following three distance measures between the normalized Fiedler vectors ⃗( 1 ) and ⃗( 2 ). First, we used the Euclidean norm between Fiedler vectors, which is given by ( ( 1 ), ( 2 )) = √∑( ( 1 ) − ( 2 )) 2 =1 . (6) Second, we used the maximum norm between the Fielder vectors, which is given by Third, we used the cosine dissimilarity between the Fiedler vectors, which is given by We used the appropriate orientation of the Fiedler vectors to calculate these distance measures (see the Supplementary Material S3 for more details).
To obtain an RP for each distance measure, one needs to define a threshold . The value of may have a crucial impact on the structure of the RP (Marwan et al., 2007). We used such that the density of points in the RPs was fixed. In other words, all RPs had the same ratio of the number of recurrences to ( − 1), the total number of possible recurrences. (Recurrence points along the main diagonal are ignored because they are trivial.) The advantage of this choice is the opportunity to compare the structure of different RPs, because such comparisons are only meaningful if the RPs have the same density of points (Marwan et al., 2007). We ran our analysis for three different thresholds such that the density of points was 0.01, 0.05, and 0.10, respectively.

Recurrence quantification analysis
To quantify the structure of the RPs, we employed 12 RQA measures (Marwan et al., 2007), i.e., the -recurrence rate ( ) and all but two measures implemented in version 5.22 of the CRP toolbox for MATLAB provided by TOCSY: http://tocsy.agnld.uni-potsdam.de. For a full description of the measures in the CRP toolbox, see Supplementary Material S5. We excluded the recurrence rate, i.e. the density of recurrence points in a RP, because we fixed this quantity to set the threshold, , to build the RPs. We also excluded the clustering ( ) because it was generally small or undefined in our RPs due to the relatively low density of points in the RPs. Among the 11 RQA measures in the CRP toolbox, four are based on the diagonal lines of the RPs, which result from recurring sequences of adjacent functional networks: the determinism ( ), the mean diagonal line length (〈 〉), the maximal diagonal line length ( max ), and the entropy of the diagonal line lengths ( ). Three RQA measures quantify the vertical lines composed of recurrent points, which represent time intervals in which the networks do not considerably change: the laminarity ( ), the trapping time ( ), and the maximal vertical line length ( max ). Another three RQA measures assess the recurrence times, i.e., the vertical distance between recurrence points in the RPs: the recurrence time of first type ( 1), the recurrence time of second type ( 2), and the recurrence time entropy ( ). The other RQA measure assesses the RP by regarding it as a network: the transitivity ( ). These RQA measures quantify different aspects of the temporal dynamics enclosed in the RPs. We used the CRP toolbox provided by TOCSY to compute these RQA measures. These 11 measures were used to assess the MEG data set and to compare pre-and post-ictal periods in the sEEG data set.
The was used to assess peri-ictal epochs in the sEEG data set. It is given by This measure counts the number of recurrence points on diagonal lines with a distance from the main diagonal. The can be considered as a generalized auto-correlation function (Marwan et al., 2007).
To assess the variation in RQA measures between pre-and post-ictal periods in the same peri-ictal epoch in the sEEG data, we computed Δ = ( post − pre )/( pre + post ), where pre is a RQA measure calculated based on the pre-ictal RP, and post is the same RQA measure calculated based on the post-ictal RP. Note that Δ varies between -1 and 1. Values close to -1 imply pre ≫ post , values close to 1 imply pre ≪ post , and values close to 0 imply pre ≈ post .

Reduction in the number of configurations
For the MEG data set, we considered four frequency bands, two functional network measures, six network distance measures, and three threshold values. Different combinations of these four factors yield different RPs. However, a majority of these RPs are redundant. Therefore, we selected representative RPs as follows. For clarity, we define a configuration as one combination of frequency band, functional network, and distance measure; we will separately consider the threshold. For example, the combination of the alpha frequency band, PLI, and the Frobenius norm is a configuration. Some of the 4 × 2 × 6 = 48 configurations may yield similar RPs and RQA results. To assess this, we studied the relations among the 48 configurations for three randomly selected healthy controls from the MEG data set, using Pearson correlation (see Supplementary Material S4 for details). This investigation yielded a few representative configurations whose RPs we used for subsequent analysis. We carried out this assessment of configurations independently for each of the three threshold values because, as mentioned, RPs are threshold-dependent and comparisons between RPs with different thresholds are not meaningful.

Statistical methods
We used the Mann-Whitney U test to assess whether the median of each RQA measure was different between the epilepsy and control groups in the MEG data set. We applied Bonferroni-Holm correction to correct p-values due to multiple comparisons across different configurations.

Dynamic MEG functional networks
To study the dynamics of functional networks inferred from MEG data, we considered signals filtered in four frequency bands, two functional connectivity measures (i.e., AEC and PLI), and six distance measures between pairs of networks. We defined a configuration as a combination of a frequency band, connectivity measure, and distance measure. We first studied the relations between the different configurations and observed that four configurations were sufficiently representative of all the configurations (see Supplementary Material S4). All of these four configurations were in the theta band. Three of them used the AEC as connectivity measure, whereas the other one used the PLI. The three representative distance measures identified together with the AEC were the Frobenius norm ( , Eq. (3)), the spectral norm ( , Eq. (5)), and the Euclidean norm between the Fiedler vectors ( , Eq. (6)). The representative distance measure identified together with the PLI was the spectral norm. We denote these four configurations by AEC+ , AEC+ , AEC+ , and PLI+ .
We used our recurrence analysis framework to compare the dynamics of resting-state MEG functional networks between people with epilepsy and healthy controls. First, we considered the AEC+ configuration. For each individual, we obtained a sequence of 506 functional networks using AEC as connectivity measure ( Fig. 1(b)). We then computed the distance between each pair of networks using as distance measure, obtaining a distance matrix ( Fig. 1(c)). Next, we identified a threshold such that 5% of the points in the distance matrix apart from the main diagonal had a distance smaller than the threshold. By thresholding the distance matrix, we obtained an RP ( Fig. 1(d)). Figures 2(a) and 2(b) show RPs from a healthy individual and an individual with epilepsy, respectively. We then used the 11 RQA measures to compare the RPs between the healthy and epilepsy groups. For example, the leftmost panel of Fig. 2(c) compares the trapping time ( ), an RQA measure, between the two groups using the AEC+ configuration. We then tested whether the median of each RQA measure was different between the two groups using the Mann-Whitney U test. We found that the recurrence time of first type ( 1) and the recurrence time of second type ( 2) were smaller in people with epilepsy than in controls in the AEC+ configuration (see the leftmost panel of Figs. 2(d) and 2(e)).
We repeated the same analysis for the other three configurations. Overall, we observed smaller 1 and 2 in people with epilepsy than in controls in the AEC+ and AEC+ configurations; 1 was higher in people with epilepsy than controls in the AEC+ configuration; was smaller in people with epilepsy than controls in the PLI+ configuration (see Fig. 2(c)). All other differences between the two groups were not statistically significant ( -values were corrected with the Bonferroni-Holm procedure). Finally, we repeated the same analysis using RPs with densities of recurrence points of 1% and 10% and found similar results. , and PLI+ . Similarly, panels (d) and (e) show box plots for the recurrence time of first type ( 1) and second type ( 2), respectively. In (c), (d), and (e), significant differences between controls and people with epilepsy are indicated by asterisks (one asterisk represents < 0.05, and two asterisks < 0.01, Mann-Whitney U test, Bonferroni-Holm corrected). We used a density of recurrence points of 0.05 to define the threshold .
Given these differences, we tested whether RQA measures could classify individuals as to whether they had epilepsy. First, we employed the 2 measure to predict whether RPs from the AEC+ configuration were obtained from people with epilepsy or from healthy controls. Note that this was the combination of RQA measure and configuration for which the p-value was the smallest among all combinations of an RQA measure and configuration when comparing the epilepsy and healthy groups. We performed the receiver operating characteristic (ROC) analysis and found an AUC of 0.76, sensitivity of 0.58, and specificity of 0.88 (see Fig. 3(a)). We then used MATLAB's Classification Learner Toolbox to classify the two groups using the 2 values in the AEC+ configuration. We found an accuracy of 69.2% using a logistic regression, i.e. 36 out of 52 individuals were correctly classified (see the blue bar in Fig. 3(b)). Next, we used the 11 RQA measures altogether and found that the classification accuracy did not improve for any of the four configurations, i.e., either 65.4 or 69.2% classification accuracy (see the black bars in Fig. 3(b)). Finally, we tested whether combining all RQA measures from the four configurations yielded a higher accuracy. Because in this case we would be attempting a classification of 52 individuals using a relatively large number of features (44 features from 4 configurations × 11 RQA measures), we first reduced the number of features by using principal component analysis. Not to bias the principal components towards RQA measures with higher variances, we normalized the features. We then observed that the first 12 principal components explained 85% of the variance of all RQA measures and all configurations. Therefore, we used them to perform the classification. These principal components yielded a similar accuracy to that for the other classification methods (see red bar in Fig. 3(b)). To avoid overfitting, we employed a 50fold cross-validation procedure in all these classifications. The best classifiers were the cosine k nearest neighbor (kNN) for the AEC+ and AEC+ configurations, the coarse tree for the AEC+ configuration and PCA, and the fine kNN for the PLI+ configuration.
We also performed the same classification analysis by (i) using the weighted mean degree of the functional networks and (ii) by applying RQA to traditional RPs obtained from the MEG time series, rather than from the dFNs (see Supplementary Material S6 and S7). We obtained a classification accuracy of 67.3% using the weighted mean degree and an accuracy of 69.2% using the traditional RQA analysis. Figure 3: Classification of people with epilepsy using RQA applied to dFNs inferred from MEG data. (a) Receiver operating characteristic (ROC) analysis to classify people as either healthy or with epilepsy using the recurrence time of second type ( 2) from the configuration AEC+ . The circle represents the optimal operating point of the ROC curve, for which the sensitivity is equal to 0.58, and the specificity is equal to 0.88. The area under the curve (AUC) is equal to 0.76. (b) Accuracy of the classification using different features from the RQA analysis and 50-fold cross-validation: the recurrence time of second type ( 2) from the AEC+ configuration (blue bar), all RQA measures from each of the four representative configurations (grey bars), and 12 principal components explaining 85% of the variance of all RQA measures across the four configurations (red bar).

Dynamic sEEG functional networks
We also applied our framework to dFNs inferred from sEEG data from people with drug-resistant focal epilepsy. In contrast to the resting-state MEG data, the sEEG data contained electrographic seizures. We considered only one broad frequency band, only the Pearson's correlation as connectivity measure, but the same six distance measures between functional networks. As a result of applying the reduction of configurations (see Supplementary Material S4), we focused our analysis on three distance measures: , , and .
We analyzed sEEG data from 10 individuals, each of them with two to four peri-ictal epochs. There were 29 peri-ictal epochs in total. Each peri-ictal epoch contained five minutes of data before seizure onset (pre-ictal), the seizure (ictal), and five minutes of data after seizure offset (post-ictal). Our aim was to observe whether functional networks show stereotypical dynamics throughout seizures, whether we can detect seizures, and how the pre-ictal, ictal, and post-ictal networks relate to each other.
For each peri-ictal epoch of each individual, we computed a sequence of functional networks using Pearson's correlation, where varied between 1542 and 1822 with median 1632, depending on the length of the ictal periods. Below, we present results obtained using the Frobenius norm, Eq. (3), to compute the distance between functional networks, and using a threshold such that the RPs had a density of recurrence points equal to 0.05. We obtained similar results using the two other distance measures and using the two other threshold values (i.e. thresholds such that the RPs had a density of recurrence points equal to 0.01 and 0.1).

RPs of multiple peri-ictal epochs
To compare different peri-ictal epochs of an individual, we first concatenated the sequences of functional networks from different peri-ictal epochs of the individual. Next, we computed the distance between each pair of networks in the concatenated sequence to obtain a distance matrix (see Fig.  1(c)). Then, by thresholding the distance matrix, we obtained an RP (see Fig. 1(d)). Figures 4(a) and 4(b) show RPs for two individuals. Each RP contains three peri-ictal epochs. We observe that most of the recurrence points are located in the pre-ictal periods (i.e., in the first, fourth, and seventh diagonal blocks) and that there is also a high density of recurrence points in off-diagonal blocks corresponding to the cross relation between different pre-ictal periods. The ictal periods also show high density of recurrence points within the same peri-ictal epoch and between different peri-ictal epochs. In contrast, the post-ictal periods show a low density of recurrence points both within and between epochs. These results imply that pre-ictal and ictal functional networks are more similar between themselves than post-ictal functional networks. There is also a considerable frequency of recurrence between pre-ictal and ictal networks, both within and between epochs. This result suggests that the functional networks involved in seizures emerge before seizure onset. In contrast, the functional networks after seizure offset are relatively different from networks during both pre-ictal and ictal periods.
To quantify the observations made from Figs. 4(a) and 4(b), we measured the relative recurrence rate ( ) defined as the fraction of the actual recurrence points in a block divided by the fraction of recurrence points expected by chance, i.e., 5% (because the overall density of recurrence points was set to 5%). Figure 4(c) shows the relative of each type of block in the RPs for all individuals. The two leftmost bars in each type of block in Fig. 4(c) correspond to the two individuals whose RPs are shown in Figs. 4(a) and 4(b) and confirm our previous observations. The relative of these two individuals further shows that, whereas the recurrence rate between a pre-ictal period and an ictal period belonging to different peri-ictal epochs is not higher than chance in Fig. 4(a), it is higher than chance in Fig. 4(b). Overall, Fig. 4(c) shows that the pre-ictal periods have recurrence rates higher than chance in 9 out of the 10 individuals. The ictal periods also have higher recurrence rates than chance in all but one individual. Post-ictal periods have low recurrence rates, except in two individuals. The cross relations between different pre-ictal periods, between pre-ictal and ictal periods, and between ictal periods show considerable variability in terms of the relative , which is either below or above chance depending on the individual. In contrast, the cross relations between pre-ictal and post-ictal periods, between ictal and post-ictal periods, and between different post-ictal periods is consistently lower than chance, except in one individual that has high relative in the cross relation between different post-ictal periods. These results support that the observations made from Figs. 4(a) and 4(b) for two individuals generalize to most of the 10 individuals. The individual represented by the third bar from the left in Fig. 4(c) (shown in dark green) is an outlier. in nine types of blocks in the RPs: the label 'pre-ictal' corresponds to the three pre-ictal diagonal blocks of the RP; 'ictal' corresponds to the three ictal diagonal blocks; 'post-ictal' corresponds to the three post-ictal diagonal blocks; 'pre|pre' corresponds to the six off-diagonal blocks that compare different pre-ictal periods; 'pre|ictal' corresponds to the 18 off-diagonal blocks that compare pre-ictal and ictal periods; 'pre|post' corresponds to the 18 off-diagonal blocks that compare pre-ictal and post-ictal periods; 'ictal|ictal' corresponds to the six off-diagonal blocks that compare different ictal periods; 'ictal|post' corresponds to the 18 off-diagonal blocks that compare ictal and post-ictal periods; and the 'post|post' corresponds to the six off-diagonal blocks that compare different post-ictal periods. For each type of block, we plot 10 bars, one for each individual. The two leftmost bars correspond to the RPs shown in (a) and (b), respectively. The chance level, i.e., 1, is represented by the solid horizontal line. Figure S2 shows that the results in (c) remain similar for other configurations (i.e. using different distance measures between networks).

RPs of single peri-ictal epochs
We then assessed whether or not RQA may be able to detect seizures. RPs comprising multiple periictal epochs are not appropriate for this purpose, because a peri-ictal epoch may have a disproportionately high or low fraction of recurrence points compared to other peri-ictal epochs for the same individual. Therefore, here we did not concatenate the sequences of functional networks over different peri-ictal epochs. Instead, we constructed an RP for each of the 29 peri-ictal epochs.
Figures 5(a) and 5(b) show RPs from a peri-ictal epoch of different individuals. These figures show that the pre-ictal periods present a higher recurrence rate than both ictal and post-ictal periods, in agreement with Fig. 4. Transitions in the RPs are also noticeable at the seizure onsets. Furthermore, whereas there is recurrence of pre-ictal networks in the ictal period, there is neither recurrence of pre-ictal nor ictal networks in the post-ictal period. Additionally, the post-ictal period has a small number of recurrence points.
To assess the potential of RQA to detect seizure onset, we computed the -recurrence rate ( ) for the RPs shown in Figs. 5(a) and 5(b). The -recurrence rate quantifies the frequency of recurrence when pairs of functional networks are time points apart. When is larger than the seizure onset time, measures recurrences of pre-ictal networks in the ictal and post-ictal periods, as well as recurrences of ictal networks in the post-ictal period. When is larger than the seizure offset time, only uses recurrences of pre-ictal networks in the post-ictal period. In contrast, when is smaller than the seizure onset time, also accounts for recurrences within the pre-ictal period and within the post-ictal period, as well as within the ictal period when is also smaller than the duration of the seizure. Because the relative is large in the pre-ictal period, but small in the cross relation between pre-and post-ictal periods, we expect that is large if and only if is smaller than the seizure onset time. Furthermore, we expect that decreases as increases from zero towards the seizure onset time because the contribution of recurrences within the pre-ictal period to decreases as increases in this range. Therefore, may be a useful tool for detecting seizure onset.
The -recurrence values for the two individuals are shown in Figs. 5(c) and 5(d). As anticipated, the decreases towards 0 as approaches the value corresponding to the seizure onset time in Fig We performed the same analysis for the other 8 individuals and found that decreases to values close to zero at seizure onset or close to seizure offset in the majority of the peri-ictal epochs considered (see Fig. 6). There are three notable exceptions (see the purple line in Fig. 6(a), the blue line in Fig. 6(d), and the green line in Fig. 6(f)), which do not show a decrease in at either seizure onset or offset. Overall, the results suggest that a decrease in to values close to zero at a certain value is indicative of a seizure onset or offset at time .  Fig. 4(a). The seizure onset and offset are indicated by the red and blue lines, respectively. (b) Same as (a) but for a peri-ictal epoch of a different individual, i.e., the first peri-ictal epoch in Fig. 4(b). (c) for three peri-ictal epochs for the individual considered in (a). (d) Same as (c) but for the RP shown in (b). In (c) and (d), the black line represents the computed from the RP in (a) and (b), respectively; the other two lines correspond to the other two peri-ictal epochs of the same individual; the dashed lines indicate time lags equal to the time of seizure onset for all peri-ictal epochs; the triangles indicate time lags equal to the time of seizure offset, where the color of the triangle matches to that of the curve.

RPs of individual pre-and post-ictal periods
We have shown that the recurrence rate within the pre-ictal period is much higher than within the post-ictal period (see Fig. 4(c)). An RP for a peri-ictal epoch, which by definition contains a pre-ictal, ictal, and post-ictal period altogether (e.g., Fig. 5), does not allow us to compare recurrence features of functional networks within the pre-ictal period to those within the post-ictal period. To address this limitation, we constructed RPs containing only a pre-ictal or a post-ictal period. We did not consider RPs only containing an ictal period because the ictal periods were typically too short for RQA. Figures  7(a) and 7(b) show an RP from a pre-ictal and post-ictal period, respectively, belonging to the same peri-ictal epoch. By construction, the two RPs had the same recurrence rate, which allowed us to compare other properties of the RPs between the pre-ictal and post-ictal periods.
Figures 7(c)-(e) show the variation of , 〈 〉, and max between all pre-and post-ictal RPs within the same peri-ictal epoch for all individuals. We observe that in most peri-ictal epochs these three measures are larger in the post-ictal than pre-ictal periods (i.e., positive variation values), meaning that post-ictal RPs have typically longer diagonal lines than pre-ictal RPs. This result suggests that post-ictal dFNs may be more deterministic than pre-ictal dFNs. For some individuals, this variation is consistently positive (i.e., the post-ictal values are larger than the pre-ictal values) for all their peri-ictal epochs (e.g., the individual represented by the right-most triangles in Figs. 7(c)-(e)). Figure S3 shows the variation of the other RQA measures. The , , , max , and values tended to be larger in the post-ictal than pre-ictal periods, whereas the opposite was the case for 1, 2, and . The measures based on vertical lines ( , and max ) indicate that dFNs are more likely to be trapped in slowly changing functional networks in the post-ictal than the pre-ictal period.

Discussion
In the present study, we have proposed to use RPs and RQA to study dFNs. The framework consists of assessing the distance between functional networks at different times and define recurrences if the distance is within a threshold. By doing so, one obtains an RP on which one can perform RQA. The RQA measures inform us about the underlying dynamics of the functional networks. We applied this framework to two data sets, (i) resting-state MEG recordings from people with JME and healthy controls, and (ii) sEEG recordings containing peri-ictal epochs of people with drug-resistant focal epilepsy undergoing pre-surgical evaluation. In the MEG data set, we found that RQA measures for dynamic resting-state functional networks differed between people with epilepsy and healthy controls. In the sEEG data set, we found that pre-ictal functional networks show high recurrence not only within pre-ictal periods but also between different pre-ictal periods, and that they also recur during ictal periods. This result implies not only that functional networks involved in seizures emerge before seizure onset, but also that they are consistent across different seizures. We also observed that RQA measures were capable of detecting seizures. Finally, we observed that post-ictal dFNs are typically more deterministic and more likely to be trapped into certain networks compared to pre-ictal dFNs.
Using the MEG data set, we found significantly smaller recurrence time in people with epilepsy compared to healthy controls (for both first ( 1) and second type ( 2) and in both the AEC+ and AEC+ configurations; see Figs. 2(d) and 2(e)). This result implies that AEC functional networks recur more quickly in people with epilepsy than in controls. This finding suggests that the space of possible functional networks may be more limited in epilepsy. We speculate that a reduced repertoire of functional networks may lead to functional brain deficits. This is in agreement with the fact that cognitive deficits are commonly observed in people with epilepsy (Holmes, 2015). Future work should test the hypothesis that reduced recurrence times predict cognitive deficits. One should also further examine such findings across different configurations (see Supplementary Material S8).
We then used RQA measures to classify people as to whether they had epilepsy and found an accuracy of 69.2% (see Fig. 3). This classification power is similar to that achieved with single restingstate scalp EEG functional networks for people with idiopathic generalized epilepsy (Schmidt et al., 2016). In our study, we considered juvenile myoclonic epilepsy, which is a different type of generalized epilepsy. Future work should test whether our findings extend to other types of both generalized and focal epilepsy. Notably, the opportunity to use resting-state data to distinguish people with epilepsy from healthy people is of great clinical value because current clinical practice depends upon the observation of epileptiform discharges which are not always apparent in a scalp EEG session (Smith, 2005). Since scalp EEG is more inexpensive and available than MEG, future work should also test whether our findings generalize to dFNs inferred from scalp EEG.
We also assessed whether the combination of our approach with others achieved a superior classification accuracy. The mean number of connections in EEG functional networks has been used to differentiate between people with idiopathic generalized epilepsy and controls (Chowdhury et al., 2016). Therefore, we computed the weighted mean degree of our MEG functional networks and used it to classify the individuals (see Supplementary Material S6). We obtained 67.3% accuracy, which was slightly lower than the classification accuracy that we attained with the RQA measures on dFNs. The accuracy did not improve by combining the weighted mean degree with the RQA measures. Additionally, we also used traditional RQA measures applied to RPs obtained from the MEG time series to classify the individuals (see Supplementary Material S7). We found an accuracy of 69.2%, i.e., equivalent to the accuracy achieved using the RQA for dFNs. Again, we found that combining the two types of RQA did not improve the accuracy of the classification. This result suggests that the two types of analysis may extract similar information. However, this does not imply that the two methods are equivalent. By inferring dFNs and studying their recurrences, we are examining the recurrences of statistical dependencies between the MEG signals. In contrast, the traditional recurrence analysis looks at the recurrence of the MEG signals themselves. Another difference between the two methods is that by computing dFNs we reduced the time resolution of the recurrence analysis, i.e., from 4 ms in the traditional RPs to 400 ms in our RPs. Future work should further explore the relation between the two methods and test whether they may complement each other on various data sets.
Using the sEEG data set, we found considerable recurrence of pre-ictal and ictal functional networks in their respective periods (see Fig. 4(c)). We further observed that ictal functional networks recurred in different seizures for most of the individuals, and that there was considerable recurrence of functional networks between different pre-ictal periods and between pre-ictal and ictal periods. Such information may be useful when evaluating whether there are one or more leading networks that sustain seizures. Non-recurrence of functional networks across different seizures may be a contraindication for epilepsy surgery because multiple seizure focus may be involved. In contrast, functional networks that recur both in pre-ictal and ictal periods across different peri-ictal epochs may support epilepsy surgery and they may be used to inform where to perform the resection Lopes et al., 2017;Lopes et al., 2018;Kini et al., 2019). Future work should test whether the combination of our framework with other recent methods that use functional networks to inform epilepsy surgery and predict surgery outcome Lopes et al., 2017;Lopes et al., 2018;Kini et al., 2019) yield superior predictions. Seizure detection is highly relevant not only for seizure management (Jory et al., 2016), but also to assist neurologists in the analysis of EEG (Adeli et al., 2003). We observed that the -recurrence rate ( ) values decreased to close to zero when approached the seizure onset or offset time in 9 out of the 10 individuals (see Figs. 5(c), 5(d), and 6).
takes advantage of the fact that functional networks frequently recur within the pre-ictal period and relatively rarely within the post-ictal period. These findings show promise for the use of our framework to detect seizures. Future work should assess whether the decrease in is specific to seizures, which one can assess by additionally measuring far from seizures (i.e., in the inter-ictal periods, which we did not consider in this study because of lack of such data). Furthermore, it should be tested whether a similar accuracy of detection is possible with scalp EEG.
We also found that the dynamics of functional networks tended to be more deterministic and more frequently trapped in certain functional networks in the post-ictal period than in the pre-ictal period (see Figs. 7 and S3). We hypothesize that individuals that show consistent differences in RQA measures between pre-ictal and post-ictal periods across seizures may be particularly suited to receiving neurostimulation treatment (Morrell, 2006). Neurostimulation can be used to modulate brain activity of people with epilepsy to avoid the emergence of seizures (Morrell, 2006). Thus, we suggest that our framework may be used for finding whether an individual presents differences in dFNs between the pre-ictal and post-ictal periods that are consistent across seizures. Such a consistency supports the use of a single stimulation strategy that may be reliable across peri-ictal epochs. Furthermore, our framework also informs us of the specific differences between pre-ictal and postictal dFNs. A personalized stimulation strategy may then be designed such that it modulates the dynamics of pre-ictal functional networks into those of the post-ictal period, whereby avoiding seizures (Dalkilic, 2017).
There is a number of confounding factors to consider when assessing our sEEG results. First and foremost, each individual had different causes of epilepsy with seizures emerging from different brain regions and each individual received different electrode implantations. Consequently, functional networks from different individuals had different numbers of nodes (i.e. channels) and covered different regions of the brain. Such differences may account for some of the variability observed among individuals. Second, even for the same individual, different peri-ictal epochs contained seizures of varying lengths, and peri-ictal epochs were at different time distances from other seizures. The distance to other seizures may be an important confounding factor when comparing pre-ictal and post-ictal periods. Third, although we used 5 min before and after a seizure as a pre-ictal and postictal period, respectively, this choice was motivated by the need of a sufficient number of functional networks for our analysis, rather than a clinically motivated decision. Different definitions of pre-ictal and post-ictal periods may yield different results.
Although we applied our framework to dFNs inferred from MEG and sEEG data, in principle the framework is applicable to any kind of time-varying network or matrix. An important requirement is to have a sufficient number of networks or matrices to reliably evaluate their dynamics and recurrences (Marwan et al., 2007). Thus, the framework may not be applicable to study dFNs derived from typical fMRI experiments due to a relatively small number of time points (Hutchison et al., 2013). For this reason, we did not explore the dynamics of functional networks during seizures in the sEEG data set; seizures were not long enough for our analysis. Additionally, although we focused on epilepsy in the present study, our framework may also be suitable to studying healthy brain function or other brain disorders.
There are other computational approaches to dFNs. Common approaches include tracking certain network measures over time (Sizemore and Bassett, 2018), using hidden Markov models (Eavani et al., 2013;Sourty et al., 2016;Vidaurre et al., 2018), and considering dynamic networks as multilayer networks (de Domenico et al., 2013;Kivelä et al., 2014;Sizemore and Bassett, 2018). Other recent approaches have used distance matrices to evaluate dFNs from fMRI (Cabral et al., 2017) and dynamic correlation matrices from scalp EEG (Rosch et al., 2018). Future work should compare these different methods to our framework to find which one better characterizes dFNs in epilepsy and other contexts, and assess whether these approaches complement each other.

Supplementary Material
Text S1: Overlap between consecutive functional networks Here we describe how we determined the overlap between consecutive segments to infer dFNs. We determined the amount of overlap as a compromise between two considerations; first, the larger the overlap, the larger is the number of functional networks that are inferred from an epoch of data, and in turn the larger and potentially richer is the resulting RP; second, an excessively large overlap produces trivial recurrences between functional networks adjacent in time that are not of interest (see section 2.4).
For a healthy control randomly chosen from the MEG data set, we first computed time-varying functional networks from the MEG signals filtered in the theta band using AEC (see section 2.3 and Supplementary Material S2). We divided the data into segments with maximal overlap. In other words, consecutive segments of 500 samples shared 499 samples (i.e., sliding the window only one sample from one segment to the next). This procedure yielded a sequence of functional networks ( ) = ( ( )), where = 1, … , and was the number of networks, i.e., the number of segments. We then employed the spectral norm distance measure, Eq. (5), to compute the distance between each pair of networks and hence a distance matrix (see Fig. S1(a)). Note that elements close to the main diagonal of the distance matrix correspond to large overlaps between pairs of functional networks, whereas elements far from the main diagonal correspond to small or zero overlaps. By definition, all main diagonal elements of this matrix are equal to zero because they correspond to pairs of identical networks. Elements close to the main diagonal were close to zero, whereas elements far from the main diagonal oscillated between different distance values. Therefore, the present RP would show a high density of recurrence points close to the main diagonal. Such recurrence points correspond to the so-called tangential motion and would not provide information about actual recurrences in the dynamical system (Marwan et al., 2007). Therefore, in RPs that we use to perform RQA, we should choose an overlap small enough to be able to avoid these meaningless recurrence points, but simultaneously large enough not to miss information about the temporal evolution of the functional networks.
To this end, we plotted the distances between pairs of functional networks against the amount of the overlap in Fig. S1(b). All the − 1 pairs of networks ( ( ), ( )) with | − | = 1 corresponded to the maximal overlap (499/500). The distance between each of these 499 pairs of networks is represented as a dot located at the overlap value of 499/500 × 100% = 99.8% in Fig. S1(b). All the − 2 pairs of networks with | − | = 2 corresponded to an overlap of 498/500 × 100% = 99.6%. In general, all pairs of networks with | − | = and < 500 corresponded to an overlap of (500 − )/500 × 100%. (For ≥ 500 the overlap was zero.) Figure S1(b) suggests that an 80% overlap between consecutive segments represents a good compromise because the minimum distance between pairs of networks, which would be relevant to recurrence events, does not substantially increase if we further decrease the overlap. The RP based on the 80% overlap is shown in Fig. S1(c). We observe that the RP avoids substantial tangential motion.
We repeated the same analysis using the other frequency bands (i.e., alpha, beta, and gamma), both the AEC and PLI as a functional connectivity measure, and the six network distance measures described in section 2.4. The results obtained for each combination were similar to those presented in Fig. S1. We also repeated this analysis for two other randomly chosen healthy controls and found similar results. Therefore, we decided to use the 80% overlap in the analyses presented in the main text. Figure S1: Assessment of the overlap between consecutive segments for the recurrence analysis. (a) Distance matrix based on the spectral norm ( ), Eq. (5). A point ( , ) in this matrix is the spectral norm of the difference between networks inferred at the time windows and . Points in between the two solid lines parallel to the main diagonal correspond to overlaps higher than 80%. (b) Spectral norm of the difference between pairs of networks as a function of the overlap. At large overlaps the distance is small, demonstrating the similarity between the networks. The dashed line at 80% overlap represents our choice. (c) Recurrence plot using the spectral norm and an overlap of 80%. Note that the temporal scale in panel (c) is different from that in panel (a); in panel (a) one time unit corresponds to 4 ms (i.e., overlap of 499/500), whereas in panel (c) one time unit is 400 ms (i.e., overlap of 400/500). We obtained these results using MEG data from one healthy individual filtered in the theta band. We employed the AEC to calculate the functional connectivity.

Text S2: Functional network measures
To construct functional networks from the MEG data, we used the phase lag index (PLI) (Stam et al., 2007) and the amplitude envelope correlation (AEC) with orthogonalized signals (Hipp et al., 2012).
Consider two sources corresponding to two nodes in the network. The PLI between the two nodes is given by where Δ ( ) is the instantaneous phase difference between the two time series at the two nodes and at time , = 1, … , ( is the total number of time points in the segment; we set = 500), sign is the sign function, and 〈•〉 is the average of the instantaneous phase difference over time.
To compute the AEC between pairs of source signals, we first orthogonalized the source signals to avoid spurious correlations due to source leakage (Hipp et al., 2012). The orthogonalization is given by where ( , ) and ( , ) are two source signals in a frequency band , * ( , ) is the complex conjugate of ( , ), ⊥ ( , ) is the orthogonalized ( , ) with respect to ( , ), and imag( ) stands for the imaginary part of . The orthogonalization is also done in the opposite direction, from ( , ) to ( , ), yielding ⊥ ( , ). We then calculated the Pearson's correlation between ⊥ ( , ) and ( , ), and between ( , ) and ⊥ ( , ), and took the average of the two values to make the AEC symmetric.

Text S3: Orientation of the Fiedler vectors
Algorithms to compute the normalized Fiedler vector ⃗( ) do not differentiate between ⃗( ) and − ⃗( ). When comparing two Fiedler vectors to calculate the distance between two networks, their orientation is not arbitrary because opposite orientations lead to different distance values. To address this, we determined the orientation of the Fiedler vectors under the assumption that if two networks are similar, then their respective Fiedler vectors are likely to be similar.
First, we computed functional networks with the maximal overlap between consecutive segments of data, such that consecutive networks shared 499 samples out of the 500 samples. Due to the high overlap, the consecutive networks were expected to be similar to each other. Second, we computed the Fiedler vector ⃗(1) of the functional network at time = 1, and chose an arbitrary orientation. Third, we computed two possible Fiedler vectors with opposite orientations for the functional network at time = 2. Fourth, we calculated the Euclidean distance between ⃗(1) and the two possible Fiedler vectors at = 2. Fifth, we selected the orientation of the Fiedler vector at time = 2 such that the Fiedler vector with that orientation was the closer to ⃗(1) in terms of the Euclidean distance. We refer to the vector with the selected orientation as ⃗(2). Sixth, we computed the two possible Fiedler vectors at time = 3 and decided its orientation by taking the vector that was the closer to ⃗(2) in the Euclidean distance. We refer to the selected Fiedler vector as ⃗(3). We repeated these steps to sequentially and uniquely determine the Fiedler vectors ⃗(2), ⃗(3), … , ⃗( ), where is the number of all functional networks. Note that our subsequent analysis only used a fraction of these Fiedler vectors. As described in the Supplementary Material text S1, we used an overlap of 80%, i.e. a shift of 100 samples out of 500 samples from one segment to the next. Therefore, we used the Fiedler vectors ⃗(1), ⃗(101), ⃗(201), and so forth.

Text S4: Reduction in the number of configurations
To remove redundant configurations from our analysis, we studied the relations between the different configurations using MEG data from three arbitrarily selected healthy controls. First, for each configuration, participant, and threshold value, we computed the RP and then the 11 RQA values. Therefore, for each participant and threshold value, we obtained a matrix of size (48 × 11) that stored the 11 RQA values for each of the 48 configurations. Then, for a given participant and threshold value, we computed the Pearson's correlation coefficient between every pair of rows of the matrix , i.e. the correlation between the two configurations by regarding the RQA values as samples. All pairs of configurations with the same functional connectivity measure and network-distance measure but different frequency bands had correlations larger than 0.9. In other words, different frequency bands yielded similar information. This result was consistent across the three participants and threshold values. Therefore, we decided to focus on only one frequency band. We chose the theta band because the original MEG recordings had the highest power in this band. By making this choice we reduced the space of configurations from 48 to 12 (i.e., 2 functional connectivity measures times 6 distance measures).
We then searched for a set of representative configurations such that all configurations had a correlation higher than 0.9 to at least one of the configurations in the set. It turned out that a set composed of at least four configurations was necessary to satisfy this condition. Finally, we searched for the set of configurations with higher correlations to the other configurations across the three controls and the three thresholds, and found a set with the following four configurations: (i) theta band + AEC + Frobenius norm ( ), (ii) theta band + AEC + spectral norm ( ), (iii) theta band + AEC + Euclidean distance between Fiedler vectors ( ), and (iv) theta band + PLI + spectral norm ( ).
In the case of the sEEG data set, we initially had 6 configurations corresponding to the 6 distance measures. Following the results for the MEG data set, we focused on 3 configurations, i.e., those based on the Frobenius norm, the spectral norm, and the Euclidean distance between Fiedler vectors.

Text S5: RQA measures
We used 11 RQA measures to quantify RPs (Marwan et al., 2007;Marwan et al., 2015). For consistency with Marwan et al.'s notation, here we use × as the size of the RP, and the indices and to denote the entries of the recurrence matrix (i.e. , instead of 1 , 2 as in the main text).
Four measures were based on diagonal lines in the RP: (1) The determinism ( ) is given by where is the length of a diagonal line in the RP (i.e., a consecutive sequence of matrix entries equal to 1 that are parallel to the main diagonal in the recurrence matrix). The quantity 1 ( ) is the number of diagonal lines of length in the RP and is given by where , is the ( , )th entry of the recurrence matrix, Eq.
(2). The is therefore the ratio of recurrence points that form diagonals in the RP of at least length min to all recurrence points. We used min = 2. Note that a diagonal line of length means that pairs of dFNs at different times remained similar to each other for the duration of consecutive time points.
(2) The mean diagonal line length (〈 〉) is defined as (3) The maximal diagonal line length ( max ) is the longest diagonal line in the RP, i.e., where { } is the set of all diagonal line lengths in the RP, and is the total number of diagonal lines, which is given by (4) The entropy of the diagonal line length ( ) is given by where ( ) = 1 ( )/ , i.e. the probability to find a diagonal line of length . Thus, is the Shannon entropy of ( ) and quantifies the complexity of the RP in terms of diagonal lines.
Additionally, we used the following three measures based on vertical lines in the RP: (1) The laminarity ( ) is analogous to the determinism and is defined as the ratio of recurrence points that form vertical lines to all recurrence points, i.e.
where is the length of vertical lines in the RP, 2 ( ) is the number of vertical lines of length in the RP, which is given by and min is the minimum length of vertical lines to be considered. We used min = 2.
(2) The trapping time ( ) is the average length of vertical lines and is given by (3) The maximal vertical line length ( max ) is defined as where { } is the set of all vertical line lengths in the RP, and is the total number of vertical lines given by (S13) Furthermore, we considered three measures that assess recurrence times. For a given column in the recurrence matrix, the elements , = 1, 1 ≤ ≤ , are its recurrence points. A recurrence time of first type of column is the number of time points from one recurrence point to the next along the column. Consecutive recurrence points, i.e. , = 1 and , +1 = 1, are counted as a recurrence time of first type equal to 1. To calculate the overall recurrence time of first type of the RP, denoted by 1, we identify all recurrence times of first type across all columns and then average them.
Recurrence times equal to 1 may result from tangential motion and not actual recurrences of the dynamical system. To account for this, one can alternatively define recurrence times as the length of time between recurrence points neglecting recurrence times of 1. The recurrence time of second type, denoted by 2, is the average of all such recurrence times in the RP. (S15) Text S6: Classification of people with epilepsy and healthy controls using the weighted mean degree We asked whether the weighted mean degree of the dFNs, which is a simpler measure of functional connectivity than those based on RQA measures, is capable of classifying people with epilepsy and controls. As we did in the RQA-based classification shown in the main text, we considered the dFNs obtained in the theta band. For each individual and each functional connectivity measure (i.e., AEC or PLI), we calculated the weighted mean degree of each of the 506 functional networks. The weighted mean degree 〈 ( )〉 at time was given by ( (S17) We then compared the values between the two groups of individuals and observed that they were not statistically different for each functional connectivity measure (i.e., AEC or PLI) (Mann-Whitney U test). Finally, using the two values as features to classify the two groups, we found an accuracy of 67.3%; we used the cosine kNN classifier which was the best classifier as identified by MATLAB's Classification Learner Toolbox.
Text S7: Classification of people with epilepsy and healthy controls using the RQA on RPs directly applied to multivariate MEG time series We compared our framework with the traditional RPs and RQA, i.e., those applied to MEG time series from individual sources instead of to dFNs. Because we only considered the theta band in the main text, we also restricted the traditional approach to the time series in the theta band. For each individual, we computed a traditional RP using Eq. (1)  (S18) where ⃗( ) = ( 1 ( ), 2 ( ), … , 90 ( )) ⊤ and ⃗( ) = ( 1 ( ), 2 ( ), … , 90 ( )) ⊤ were vectors from the source reconstructed signals at times and , respectively. Note that, whereas in the dFNs framework we used 500 samples to generate a network and an overlap of 80% between consecutive networks, here we used every sample as an independent time point. Therefore, the traditional RP was much larger than an RP for dFNs. For computational tractability, we restricted the size of the RP to 5000 × 5000 constructed from the first 5000 samples. We set the threshold on the distance values to define recurrence to 0.05, which was the same value as the one used in the main text. Next, we used the RQA measures to quantify the RPs of each individual. Finally, we employed the RQA measures to classify people as healthy or having epilepsy and obtained an accuracy of 69.2%; we used the linear support vector machine which was the best classifier as identified by MATLAB's Classification Learner Toolbox.
Text S8: A note on the RQA measures that differ between people with epilepsy and healthy controls in the MEG data set We found that both the recurrence times of first ( 1) and second type ( 2) are smaller in people with epilepsy compared to healthy controls. However, this finding was specific to the AEC+ and AEC+ configurations. In the AEC+ configuration we found a significantly higher 1 in people with epilepsy but no statistical difference in 2 between the two groups. Given that the difference between 1 and 2 is that 1 is affected by tangential motion, this finding implies that RPs from controls in the AEC+ configuration reflect more tangential motion than those from people with epilepsy. We also found that the trapping time ( ) was significantly smaller in people with epilepsy than controls in the PLI+ configuration (see Fig. 2(c)). The is the average length of vertical lines in the RP and measures the time during which dFNs are trapped near a certain network. In this way, PLI networks from healthy controls are more likely to be trapped than those from people with epilepsy. In contrast, the of AEC networks is not significantly different between the two groups in any of the three considered configurations. The fact that we observed different results for different configurations is not a contradiction because different configurations assess different features of the functional networks. Figure S2: Relative recurrence rate ( ) for the different blocks of the RP from sEEG data. Panel (a) and (b) are equivalent to Figure 5, except that to compute the relative , we constructed the RPs using (a) the spectral norm, Eq. (5), and (b) the Euclidean norm between Fiedler vectors, Eq. (6). The 9 different types of blocks in the RP are described in the caption of Figure 4(c). As in Figure 5, for each type of block, we plot 10 bars of different colors, one for each individual. To obtain the RPs, we used a density of recurrence points of 0.05. Figure S3: Variation in RQA values between pre-and post-ictal RPs of sEEG dFNs. Each panel represents the variation in one RQA measure, as in Figs. 7(c)-(e). In each panel, each triangle corresponds to a different peri-ictal epoch, and the triangles of the same color belong to the same individual. To compute the underlying RPs, we used the Frobenius norm as a pairwise distance between networks, Eq. (3), and a density of recurrence points of 0.05 to define the threshold .