Connectivity detection in application to spike-wave discharge study

In our study, we compare three popular approaches to directed coupling analysis, in particular transfer entropy and two types of Granger causality, applied to real data from genetic absence epilepsy rats. We have chosen the channels for which the coupling architecture is already well known from previous studies. Recordings from 5 WAG/Rij rats of 8 hours duration with at least 28 spontaneous seizures of length not less than 6 s in each recording were studied. To test results for significance, surrogate signals based on series permutation technique were constructed. Connectivity development in time was investigated by considering six two-second intervals before, during and after the seizure. Our outcomes showed large differences between studied approaches, while all of them exploit the same general idea. Transfer entropy demonstrated the smallest number of significant couplings throughout all three considered measures, while the linear Granger causality showed the largest number of them. This indicates that transfer entropy is the most conservative measure and the least sensitive one. Its sensitivity is affected by insufficient series length. The linear Granger causality is likely to demonstrate insufficient specificity.


Introduction
In this paper, we aim to investigate the variability of coupling estimation results using different directed approaches in application to brain connectivity study. Most papers, where new coupling measures were proposed, also provided some comparison with previous studies. The applications to EEG data were performed in many such papers, see e. g. [Baccala and Sameshima, 2001;Sysoeva et al., 2014;Sommerlade et al., 2015]. In addition, there is a number of reviews on coupling method comparison till now, including those which consider primary neurophysiological applications, see [Pereda et al., 2005;Gourévitch et al., 2006;Bezruchko et al., 2008]. These comparisons were mostly done based on simu-lated data, since such an approach provides full control, all couplings are reliably known and may be set by hand. The idea of such an approach is quite clear: first, to study the method reliability on simulated examples, and second, to get some new results from experimental data. Here, we propose a different approach: to use well studied physiological system for which the couplings are mostly known based on both anatomy and physiology. We consider cortico-thalamo-cortical network involved into absence seizures as such a system.
Absence epilepsy is a widespread form of epilepsy, providing up to 50% of cases among children and adolescents [Megiddo et al., 2016]. The etiology of this form is not completely clear, being classified as genetic by ILAE (International league against epilepsy). The main manifestation of this epilepsy is partial or complete loss of consciousness for short time intervals [Coenen and van Luijtelaar, 2003;Volnova and Lenkov, 2012] -the duration of the seizure is, as a rule, about 5-10 s [Holmes et al., 1987;Bosnyakova et al., 2007;Akman et al., 2010]. Spike-wave discharges (SWDs) are the main manifestation of absence epilepsy on the electroencephalogram (EEG), including both scalp EEGs and recordings of local field potentials (LFPs) from cortex and thalamus.
Due to the fact that children do not have clinical indications for intracranial surgery, and due to the fact that the skull works as a filter for brain signals surface electroencephalograms (EEG) do not carry enough information about the activity of the involved brain structures. This applies primarily to the thalamus, which is responsible for the epileptic seizure generation according to the modern concepts [Meeren et al., 2005]. Therefore, the main results in the investigation of absence epilepsy were obtained on genetic models -mainly rats, because access to the deep brain structures (primarily the thalamus) is necessary to study mechanisms of generation and maintenance of absences [Russo et al., 2016]. The works [Meeren et al., 2002;Lüttjohann and van Luijtelaar, 2012;Sysoeva et al., 2016a;van Rijn et al., 2010] used a very large data sets: at least 16 animals in each experiment and several hundred seizures in total (usually not less than 10 seizures per animal). Different methods were applied to study the problem, including linear and nonlinear correlation analysis [Meeren et al., 2002;Lüttjohann and van Luijtelaar, 2012], mutual information function [Grishchenko et al., 2017], synchronization measures [Maksimenko et al., 2017], linear [Sitnikova et al., 2008] and nonlinear Granger causality [Sysoeva et al., 2014;Sysoeva et al., 2016a]. So, the connectivity in cortico-thalamo-cortical circuit at absence epilepsy seems to be very well understood.
Since we aim to track coupling changes over time to follow the processes in coupling leading to development of SWDs, the short time series have to be used. This actually means that lack of data prevents the use of multivariate approaches, such as conditioned Granger causality of partial directed coherence [Baccala and Sameshima, 2001], as well as methods based on phase dynamics modeling [Rosenblum and Pikovsky, 2001;Smirnov and Andrzejak, 2005;Navrotskaya et al., 2019] which usually demand at least 30-50 characteristic periods in order to obtain reliable results [Smirnov and Bezruchko, 2003]. Therefore, here we focused on pairwise approaches resting the consideration of multichannel techniques for future.

Experimental Data
All experiments were approved by the Radboud Ethics Committee of the University of Nijmegen (RUDEC 2006-064). Animals were operated on under isoflurane anesthesia. The electrodes were installed in accordance with the atlas [Paxinos and Watson, 2006] with the following coordinates: frontal cortex (FC) [AP +3.5; L3], parietal cortex (PC) [AP −1.6; L4], and occipital cortex (OC) [AP −6; L −3.5] and inserted in the hippocampus (HP) [AP −3.5; L2; depth: 3.5]. The reference electrode was placed in the cerebellum. Histological investigation of electrode position was performed and only animals with verified and correct electrode positions were included in the data analyses. After the electrodes were installed, the animals were placed in separate cuvettes, they received enough water and food, and were rehabilitated for at least 14 days before the data were collected. These channels were chosen following the goal of our study. When the study of efficiency of different techniques is performed on simulated data, usually unidirectionally coupled, bidirectionally coupled and uncoupled channel pairs have to be considered to answer the questions: 1) whether the studied technique is able to detect the actual coupling, 2) whether it considers the absent coupling as insignificant or weak, 3) if it is able to distinguish between unidirectional and bidirectional coupling. In real world applications we cannot be completely sure, but based on known results [Sysoeva et al., 2016b] hippocampus is not involved into SWDs, PC and FC are coupled bidirectionally and OC is mostly passively involved in seizures, obtaining inputs from PC and FC.
The 130-minute recordings of intracranial EEG (signals of local field potentials, LFPs) from 5 male WAG/Rij rats were collected. All animals had at least 28 spontaneous seizures of length not lesser than 6 s in each recording, with all electroencephalographic channels demonstrating fine signal amplitude along whole recording.
The data were digitized with sampling frequency of 512 Hz and recorded with 16-bit analog-digital converter with hardware-based filtering in the range 1-99 Hz, including suppression of 50 Hz. The same data was previously partially used in a pharmacological experiment [van Rijn et al., 2010] as a control group. The analysis of mutual information between channels in a moving window was carried out in the work [Sysoeva et al., 2016a;Sysoeva et al., 2016b], but only animal-averaged curves were analyzed. Each animal was considered individually. The length of the selected seizures was at least 6 s, and it was taken five seconds before and after the seizure (no SWDs occurred in these 5 s intervals). From each seizure, six 2-seconds time intervals were selected (including baseline, preictal and postictal activity epochs), and the coupling measures were calculated for which epoch separately. For convenience, the seizure onset was taken as 0 and the seizure termination was taken as T . Then, the studied intervals, selected in accordance with the results of the work [Sysoeva et al., 2016a], can be designated as follows (see Fig. 1 for visualization): 1. [−5; −3] s before the seizure onset (background activity); 2. [−2; 0] s before the seizure onset (preictal activity); at this time, as a rule, changes in connectivity leading to seizure initiation are detected already [Sysoeva et al., 2014;Lüttjohann and van Luijtelaar, 2012]; 3. The time series of recording intracranial EEG for the parietal cortex is shown in Fig. 1 and the spectra for all 6 segments averaged over all 28 seizures are depicted in Fig. 2. In Fig. 1, interval [0; 9] s after the seizure onset has a larger amplitude compared to the background [−5; 0] and [9; 12] s. This behavior is typical for absence epilepsy [Lüttjohann and van Luijtelaar, 2015].

Granger Causality
The idea of the Granger causality method is as follows. Let there be two objects: an object X, from which the time series is derived {x n } N n=1 , and an object Y , from which the time series is derived {y n } N n=1 . Studying the causal interactions between X and Y using Granger causality involves three steps. First, a univariate predictive model is constructed based on the data from only one series {x n } N n=1 , the effect on which is estimated: where x n is the predicted value corresponding to the measured value x n ; f is an approximating function, in our case it is a polynomial with order P from D s variables [Chen et al., 2004]; x n = x n , x n−l , . . . , x n−(Ds−1)l is a state vector as defined by means of the method of delays [Packard et al., 1980;Kugiumtzis, 1996]. The method of delays is a classical approach to expand time series in phase space, i. e. to obtain the high-dimensional state vector { x n } N −(Ds−1)l n=1 from the scalar time series {x n } N n=1 by shifting the series back in time (D s − 1) times with the delay (lag) l. N is the efficient length of time series {x n } N n=1 , calculated as N = N − ((D s − 1)l + τ ); τ is the length of prediction interval (prediction length), i. e. the time lag between the last point used for vector reconstruction and the predicted point. Model coefficients are estimated using the least squares method [Silverman, 1986], i. e., minimizing the square error of the prediction ε 2 s , which measures the difference between the predicted and observed values x n and x n correspondingly.
where σ 2 s is the empirical variance of the time series {x n } N n=1 . At the second stage, the bivariate model (3) is built on the basis of both time series {x n } N n=1 and {y n } N n=1 : x n+τ = g x n , x n−l , . . . , x n−(Ds−1)l , y n , . . . , y n−(Da−1)l , where D a is the dimension of the additional part of the state vector y n = y n , y n−l , . . . , y n−(Da−1)l constructed from a scalar time series {y n } N n=1 , and x n is the value predicted with the bivariate model and corresponding to x n . Thus, the dimension of the bivariate model can be calculated as D j = D s + D a , and the prediction error is denoted as ε 2 j and calculated analogously to (2). At the third stage, the P I prediction improvement value is calculated using (4), considered as the main characteristic of the Granger causality The situation when ε s = ε j suggests that the data from the time series {y n } N n=1 do not improve prediction of {x n } N n=1 . In other words, Y does not drive X. The situation, when ε s > 0 and ε j → 0, P I → 1, suggests that the data of the second time series {y n } N n=1 significantly improves prediction of the first one {x n } N n=1 , leading to the outcome that Y drives X.
The results of the Granger causality depend on the parameters of the model [Papana et al., 2013], for example, on the type of basis function [Kornilov et al., 2016;Marinazzo et al., 2006], the polynomial order for polynomial functions [Chen et al., 2004], lag l used for state vector reconstruction and prediction length τ . Ultimately, the choice of approximating nonlinear functions in the model and its parameters are crucial for the quality of prediction and practical application of the Granger causality method, this should be considered to achieve reliable results. Univariate model, bivariate model and their parameters should be finely tuned to avoid misleading results.
When testing for causality using Granger's approach, standard linear type models (5) are most often used, as proposed by [Granger, 1969] and considered by [Gourévitch et al., 2006].
where c s i and c j i are coefficients for univariate and bivariate models correspondingly usually fitted with leastsquares routine.
Nonlinear models can provide significantly better sensitivity, as it was multiply shown [Chen et al., 2004;Marinazzo et al., 2006;Sysoev and Sysoeva, 2015]. However, its specificity could be insufficient, providing a lot of false positive outcomes. Therefore, we used the model specially designed in  using nonuniform embedding to reduce the number of coefficients c s i and c j i and Bayesian information criterion by [Schwarz, 1978]. into account the value of the experimental data delayed from the predicted time point with a characteristic period T (since there is a well pronounced main time scale for absence seizures, ). The value of τ was selected following , the value of l -following [Grishchenko et al., 2020] and the value of P and D s -following [Kornilov et al., 2016]. A typical view of the dependence of the prediction improvement on time constructed for the six considered time intervals is shown in Fig. 3(a).

Transfer Entropy
The transfer entropy proposed by [Schreiber, 2000] is used to determine the directional connectivity of oscillatory systems by their time series. The idea is similar to Granger causality; but instead of construction of predictive models like (5) and (6) one estimates conditional distribution functions. It was shown by [Barnett et al., 2009] that for some very simple cases transfer entropy and Granger causality can be equivalent. However, in most cases of interest such equivalence cannot be shown or it is not the case.
Estimating transfer entropy means actually estimating multidimensional probability densities. Therefore, many approaches designed for calculation of mutual information function can be used to calculate transfer entropy. The direct "naive" approach based on splitting the space into polygons or cubes demands a lot of data and has pure properties for nonuniform distributions providing large errors. Therefore a number of advanced methods were proposed by [Kraskov et al., 2004;Silverman, 1986;Moddemeijer, 1989;Lee et al., 2012;Darbellay and I., 1999;Kugiumtzis, 2012;Jizba et al., 2012].
The method of calculating the transfer entropy used in this work was described in the article [Sysoev, 2016] and is based on the nearest neighbors approach for calculation of mutual information by [Kraskov et al., 2004], as one of the simplest and efficient for relatively small samples (it is quite unbiased for N ≥ 10 3 ).
T E Y →X = ψ(K) + ψ(ν xn + 1) − ψ(ν xn,xn+τ + 1) − ψ(ν xn,yn + 1) n=1,N −τ , where X, y, n, N and τ are of the same meaning as in (5) and (6), K is the number of the nearest neighbor used for calculation (we used K = 6 due to short series did not allow us to use larger values as recommended in [Vlachos and Kugiumtzis, 2010;Faes et al., 2015]), ψ is a digamma function. To explain the meaning of ν let us consider (n) -the distance in the 3-dimensional state (X n , X n+τ , Y N ), which is calculated following eq. 8.
So, ν xn is number of elements in series {x n } N n=1 , which distance to a point x n strictly less than (n)/2, ν xn,yn is number of points from two-dimensional space (X n , Y n ), similar to ν xn and ν xn,xn+τ is number of points from the two-dimensional space (X n , X n+τ )) being in the neighborhood of the point (x n , x n+τ )).
The estimation of the transfer entropy by the experimental series using τ = 1, as proposed in [Kraskov et al., 2004], did not reveal the directional relationship, therefore the value τ = 6 was used similar to the adapted nonlinear approach with τ = 1.

Statistical Analysis
The estimates of all three considered measures (TE, linear and nonlinear Granger causality) varied largely for different episodes, as it is shown in Fig. 3(a) for P I calculated using nonlinear model (5). Therefore, we performed the statistical analysis to be able to compare different measures and different epochs.
To test the results of coupling analysis for significance, surrogate time series for each animal were built separately by realization permutation since the classical approach which uses phase randomization [Theiler et al., 1992] occurred to be suboptimal for data of considered type as it was shown based on simulated seizures in [Sysoev and Sysoeva, 2015]. There were 27 · 28 = 756 surrogate pairs, where 28 is the number of seizures considered for each animal. It is combinations of all possible episodes for a channel pairs, except when they are from the same episode. This number of surrogates allows one to get a high confidence level of 99.87% (p-value is equal to 1/(27 · 28 + 1) ≈ 0.0013), what is important because there is a lot of multiple testing. Further, the number of significant findings for each animal was estimated over all 28 intervals considered, see Fig. 3(b), and corrected for multiple testing.
Since we consider 5 rats, 6 time intervals, 12 channel pairs and 28 seizures for each animal, we have to make a Bonferroni or another correction to be sure that significance of results was not corrupted due to multiple testing. The number of significant results k per animal per channel pair per episode is distributed following binomial law F L,p (k) with parameters L = 28 (total number of seizures) and p ≈ 0.0013. If we set the final (after multiple testing correction) desired confidence level to be p ≤ 0.001, four or more significant values of any measure (TE, LGC or NGC) are enough to treat the result significant due to high confidence level p.
To be able to compare different number of significant results for equivalence, e. g. 14 in one case, and 22 in another, the Mann-Whitney U-test was be used. This is possible since the number of investigated cases is 28, which usually is considered to be enough.
We use commonly acceptable ideas of method sensitivity and specificity when comparing results of different techniques. We assume that the specificity is an avoidance of false positive results, i. e. ability not to detect the absent link as valid (the highest specificity means no false positives), and the sensitivity is an ability to detect really existing couplings as frequently as it is possible (the highest sensitivity corresponds to detection in 100% cases).
3 Results Fig. 5 shows that in the background (the first points on the plot at −5 s) all considered measures cannot detect statistically significant couplings. This behavior is also observed for coupling 2 seconds before the seizure onset (the second points on the plot at −2 s), and for coupling in the final stage at the seizure termination (the last points on the plot at 6 s), except that TE indicates some couplings still present in F C ↔ P C and OC → F C pairs.

Cortico-hippocampal Interactions
All three measures indicated that none of the animals had significant couplings between the hippocampus and frontal cortex during seizures. Only a single rat (No. 5) had significant connections in the pair HP→PC during normal activity (before and after the seizure). Two animals (No. 3 and No. 5) showed the coupling in the pair Hp→OC ictally. These connections were found using linear and nonlinear Granger causalities; transfer entropy did not indicate them as significant. All these connections were not very strong, see Fig. 4). The hippocampus is a part of the limbic system, which is traditionally considered to be not involved in absence seizures [Lüttjohann and van Luijtelaar, 2015]. The known evidence of its partial involvement is in another animal model of absences [Arcaro et al., 2016]. The occipital cortex is strongly anatomically connected to entorhinal cortex, another part of the limbic system. Therefore, the coupling detected between Hp and OC is not surprising. Thus, we can conclude that all methods agree that the hippocampus is practically not involved in the seizure, as it was suspected initially.

Intracortical Interactions
In Fig. 5, the calculation results of all three coupling measures are shown for a single animal: red color indicates the transfer entropy, blue color shows the linear Granger causality, and gray color shows the nonlinear Granger causality. This animal was plotted since it demonstrated more significant couplings than others. The studied intervals are shown along the X-axis. The total number of seizures for which the conclusion about coupling at a considered interval was significant is shown on the Y-axis on the left. The same number, but as a percentage of the total number of seizures (28 for all considered cases) is shown on the right. For all rats, the channel pairs for which a significant coupling is detected at the appropriate time interval are summarized in the table 1. The results are based on the statistical corrections described above. Numbers of all considered recordings are placed in the table horizontally. Six studied time intervals are placed vertically. The results are given for all three connectivity measures sequentially.
All three measures showed an increase in coupling during the ictal epochs (time points 0, 2 and 4) in the channel pairs FC↔PC and FC→OC. This effect differed in different animals (see table1); sometimes the linear and nonlinear GC demonstrated the same results, sometimes the linear GC is closer to the TE rather than to the nonlinear GC. Also, different measures showed an increase at different time points. PC→OC coupling was detected only using Granger causality, while the TE showed no significant results.
The nonlinear GC and TE did not detect significant coupling when the effects of the occipital cortex on the parietal and frontal cortex were investigated, while linear GC detected some coupling. Theoretically, OC is involved into the seizures only passively [Lüttjohann and van Luijtelaar, 2015]. So, we can conclude that the linear GC was likely to show the results in a wrong direction (not enough specificity), as it was shown in [Sysoev and Sysoeva, 2015] for linear GC.
At the seizure onset (marked 0 s on the plot and in the table 1), a statistically significant coupling is observed from the frontal and parietal cortex, and this relationship is often the most significant one (diagnosed for the largest number of seizures) all over the studied intervals. Similar behavior was observed for all investigated animals (see table 1 for details).
Generally, it can be noted that the transfer entropy showed the least number of significant coupling in comparison to two other measures, and the linear Granger causality showed the largest number. Fig. 5 shows that the transfer entropy is not sufficiently sensitive, although in some cases it also reveals a coupling.
We analyzed all animals including all channel pairs. The complete analysis of all channel pairs was summarized in the table 1, where all detected significant connections were listed for all six epochs. However, let us consider first the connections most relevant for absence epilepsy, which are between FC and PC. Common was

Discussion and Conclusion
First of all, let us note that all three considered measures (as well as some measures not considered here, like partial directed coherence [Baccala and Sameshima, 2001; Baccala et al., 2007]) are based on the same idea originating from Wiener's paper [Wiener, 1956]. Since these methods do not consider enough a priori information about investigated systems [Bezruchko and Smirnov, 2010], they can measure some effects, but would never be absolutely robust and precise, being unable to reveal mechanisms of coupling [Bezruchko and Smirnov, 2010;Barrett and Barnett, 2013]. Therefore, all outcomes in method comparison should be considered, with keeping properties of the considered data (series length, spectral characteristics, sampling frequency and so on) in the mind. Another data properties could lead to different in some kind results like outcomes of [Pereda et al., 2005]. Also, method development during last years including new approaches to coupling measure calculation provided in [Kraskov et al., 2004] and in [Marinazzo et al., 2008], and new techniques to fitting methods to data [Kornilov et al., 2016] depicted some outcomes of previous studies as not up to date. For instance, by [Gourévitch et al., 2006] the linear Granger causality was preferred over the nonlinear one due to authors were not able to interpret the results of nonlinear causality properly. However, the main source of this outcome could be inappropriate parametrization as it was shown by [Papana et al., 2013;Kornilov et al., 2016]. Also, the results of linear Granger causality application to LFP data by [Sitnikova et al., 2008] were reconsidered by [Sysoeva et al., 2014], showing different effects,

Common in Three Investigated Techniques
The common for three considered measures is as follows. First, all three measures showed very small, usually insignificant number of interactions during baseline, preictal and in postictal epochs for all rats between considered channels. This is interesting because there is al-ways some connectivity in the brain. And this means that week couplings which are typical for the normal brain cannot be reliably detected from time series consisting of 16-18 oscillations (2 s) and ∼ 10 3 data points. We know that during the seizure coupling rises a lot and only this makes methods to detect it.
Second, the methods did not detect coupling from cortical channels to the hippocampus and mostly in the opposite direction. This means that a network node which does not participate in the considered activity at all, as hippocampus in our case, is correctly detected as isolated.
Third, all free methods showed mutual increase in coupling between PC and FC. Also, all of them showed that OC is driven either by FC or by PC (sometimes by both of them) during the seizure.

Difference in Three Investigated Techniques
Nonlinear measures showed OC to be involved into seizures only passively, driven by either PC or FC, that matches the modern views on absence epilepsy. But the linear GC showed bidirectional coupling. This indicates that the linear GC method is not specific enough, being unable to distinguish between unidirectional or bidirectional connectivity as it was shown before for simulated series in [Sysoev and Sysoeva, 2015].
For each animal, the transfer entropy revealed the least number of significant links, and for the rat No. 2 -does not reveal any link at all. It can be assumed that this measure has insufficient sensitivity for the used time series length. This can be explained by the fact that the method is mainly non-parametric. It relies on estimating the distribution density in a multidimensional space, which always requires large amount of data [Schreiber, 2000;Kraskov et al., 2004].
The linear Granger causality method revealed the maximum number of significant couplings. This can be explained due to its lack of specificity, as it was shown earlier [Sysoev and Sysoeva, 2015]. Notably, this method is the only which diagnoses some couplings even in the background period.

Recommendations
The recommendations for usage of various coupling measures are as follows. For the preliminary analysis of large amounts of data (one hundred of seizures per patient/animal or more with large sampling frequency, at least 2000 Hz), one can use the transfer entropy. This is a non-parametric method that does not require painstaking selection of parameters, so it can be implemented faster than parametric ones. If a large amount of data is used, averaging should reveal the most significant interactions. To identify more subtle effects, it is desirable to use a nonlinear Granger causality method with carefully selected parameters. The linear Granger causality method, unfortunately, is difficult to recommend because of its poor specificity.