[1] 7w DYNAMICS OF DIRECTIONAL COUPLING UNDERLYING SPIKE-WAVE DISCHARGES

[1] 343w Spike and wave discharges (SWDs), generated within cortico-thalamo-cortical networks, are the electroencephalographic biomarker of absence epilepsy. The current work aims to identify mechanisms of SWD initiation, maintenance and termination by the analyses of dynamics and directionality of mutual interactions between the neocortex and various functionally different thalamic nuclei. Methods: Local-field potential recordings of 16 male Wistar Albino Glaxo from Rijswijk (WAG/Rij) rats, equipped with electrodes targeting layer 4-6 of the somatosensory cortex, rostral and caudal reticular thalamic nuclei (rRTN and cRTN), ventro-posteromedial (VPM), anterior (ATN) and posterior (PO) thalamic nuclei, were obtained. 3 s epochs prior to SWD onset, after SWD onset, prior to SWD offset and after SWD offset were analyzed with newly developed time-variant adapted nonlinear Granger causality. Results: A gradual increase in coupling toward SWD onset between cortico-cortical pairs appears as early as 2 s preictally. Next first unidirectional increase in coupling is noticed in a restricted number of cortico-thalamic and thalamocortical channel pairs, which turn into bidirectional coupling approaching SWD onset, and a gradual increase of intrathalamic coupling. Seizure onset is characterized by a coupling decrease for more than a second in a majority of channel pairs, only the cortex kept driving the cRTN. Intrathalamically the cRTN drives the PO, VPM and ATN. Most channel pairs no longer show differences in coupling with baseline during SWD maintenance, a major exception is the unidirectional coupling between cortex and cRTN. Toward the end of SWDs, more and more channel pairs show an increase in often bidirectional coupling, this increase suddenly vanishes at SWD offset. Conclusion: The initiation of SWD is due to a gradual increase in intracortical coupling, followed by a selective increase in first unidirectional and later bidirectional coupling between the cortex and thalamus and also intrathalamically. Once the network is oscillating, coupling decreases in most of the channel pairs, although the cortex keeps its influence on the cRTN. The SWD is dampened by a gradual increase in coupling strength and in the number of channel pairs that influence each other; the latter might represent an endogenous brake of SWDs.

[1] 148w Absence epilepsy is classically considered as a nonconvulsive generalized epilepsy (classification of the International League Against Epilepsy, ILAE; Berg et al., 2010) of unknown etiology. Its clinical symptoms are rather small and might even go unnoticed especially when they last only several seconds, facial automatisms might be lacking and the brief decrease of responsiveness, the interruption of ongoing behavior, and the impaired mental functioning are difficult to notice. In contrast, the electroencephalogram (EEG) during typical absence seizures is an archetypical clear and easily recognized 3-4 Hz pattern of spike-and-wave discharges (SWDs). The SWDs are conceptualized to originate at some point within the cortico-thalamic system, and rapidly engage other parts of a bilaterally distributed circuit. Sites of origin can be visualized with imaging and signal analysis techniques and nowadays most often cortical origin sites are reported (Holmes et al., 2004;Westmijse et al., 2009;Tenney et al., 2013;Ossenblok et al., 2013).

[2] 364w Similar EEG paroxysms, SWDs, appear in rat strains with a genetic predisposition to develop absence epilepsy, such as GAERS (Genetic Absence Epilepsy Rats from Strasbourg -(Vergnes et al., 1987)) and WAG/Rij (Wistar Albino Glaxo from Rijswijk - (Sitnikova and van Luijtelaar, 2006)). The EEG waveform and duration (1-30 s, mean 5-6 s) of SWD in rats and in humans are comparable, but the frequency of SWD in rats is higher 7-11 Hz (van Luijtelaar and Coenen, 1986;Sitnikova and van Luijtelaar, 2007). In WAG/Rij rats a consistent initiation zone in the perioral region of the http://dx.doi.org/10.1016/j.neuroscience.2015.11.044 0306-4522/Ó 2015 IBRO. Published by Elsevier Ltd. All rights reserved. somatosensory cortex was found (Meeren et al., 2002). This leading role of a cortical initiation zone in spreading of SWD in humans and in WAG/Rij rats was established with signal analytical methods that did take into account simultaneous processes in different locations (network analyses). Time resolved methods of coupling analysis applied for multichannel time series were used such as nonlinear association analyses, but also linear Granger methods (Granger, 1969) and phase synchronisation (Pijn et al., 1989;Westmijse et al., 2009;Sitnikova et al., 2006;Lu¨ttjohann and van Luijtelaar, 2012;Lu¨ttjohann et al., 2013). While earlier studies focused on cortical spreading and confirmed a leading role of the somatosensory cortex in the initiation of SWD (Meeren et al., 2002), subsequent studies found strong and reciprocal coupling between the frontal cortex and thalamus or increase in SWD-related increase of intra-and interhemispheric and intrathalamic coherence during SWD (Sitnikova et al. 2006(Sitnikova et al. , 2008)); the role of subparts of the thalamus was less well investigated. Classically the ventro-posteromedial (VPM) and rostral reticular (rRTN) thalamic nuclei were thought to be the primary nuclei in SWD occurrence and maintenance. The thalamus is a collection of functionally heterogeneous nuclei (sensory, limbic, motor, arousal) reciprocally connected with different parts of the cortex. Recent studies toward interactions of the assumed cortical site of origin of the SWD with various thalamic nuclei with time resolved methods revealed an additional role for the posterior nucleus (PO) of the thalamus, and for the caudal part of the reticular nucleus (cRTN) in SWD initiation and maintenance (Lu¨ttjohann and van Luijtelaar, 2012;Lu¨ttjohann et al., 2013).

[3] 106w In comparison with traditional methods of network analysis, such as cross-correlation, coherence, phase synchronisation, Granger causality (GC) may detect weak or hidden coupling, which not necessarily lead to synchronisation. It defines next to changes in coupling strength also changes in the direction of coupling within a network. Granger used only linear predictive (autoregressive) models; new nonlinear models were successfully developed and applied (Wang, 2007;Bezruchko and Smirnov, 2010). This is not trivial for its application in EEG paroxysms considering that seizure activity has nonlinear properties (Le van Quyen et al., 1999;Lehnertz, 1999;Lopes da Silva et al., 2003) and linear methods may capture only part of the coupling.

[4] 192w Here a recently developed new nonlinear approach called time-variant adapted GC was used and applied to a previously published data set (Lu¨ttjohann and van Luijtelaar, 2012;Lu¨ttjohann et al., 2013Lu¨ttjohann et al., , 2014) ) of in vivo local field potentials (LPFs) data recorded by means of intracranial electrodes implanted in the deep somatosensory cortex and in five different parts of the thalamus in WAG/Rij rats. Recent data in humans and in these genetic absence epileptic rats suggest that SWD do not arise out of the blue (Holmes et al., 2004) but are preceded by precursor and network activity in and between cortex and thalamus (Gupta et al., 2011). Also the end of SWD was initiated by a decrease of linear coupling from the somatosensory cortex to the rRTN, as well as increased coupling from the caudal to the rRTN (Lu¨ttjohann et al., 2014). Here the dynamics of corticocortical, cortico-thalamic and thalamo-thalamic network interactions will be investigated in this absence model in three stages using this new nonlinear method (Sysoeva et al., 2014): the transition from preictal to ictal phase, the ictal phase and the transition from the ictal to the postictal state.

[1] 144w Dependences of PI 0 (t) for all possible channel pairs are plotted in Figs. 2345, where Fig. 2 considers the intracortical interactions, Fig. 3 -the cortico-thalamic interaction (in direction from cortex to thalamus), Fig. 4 -thalamo-cortical interaction (from thalamus to cortex), and Fig. 5 -intrathalamic interactions. All clusters of three points significantly higher than baseline are plotted in black, others -in light gray. Red points indicate a significant decrease of PI 0 (t) after the maximum. In all cases only if there are at least three significant values in a row, they are considered as truly significant in order to avoid false positive due to repetitive testing (Maris and Oostenveld, 2007;Lu¨ttjohann and van Luijtelaar, 2012). If there are three or more points significantly different from the baseline level, this difference was established to occur in the time moment, corresponding to the first of them.

[2] 75w The preictal changes were established in all of the combinations of intracortical channel pairs characterized by a gradual increase in PI 0 (t) with a maximum immediately preceding SWD onset. The data are presented in Fig. 2 and Table 2. Earliest significant increase in PI 0 (t) was seen at or more than 2 s before seizure onset in all pairs except ctx5 ? ctx4, in which the increase occurred immediately prior to SWD onset.

[3] 30w A significant decrease in PI 0 (t) following a significant preictal increase was found for the following cortico-cortical channel pairs ctx4 M ctx5 (bidirectional), ctx6 ? ctx4, ctx6 ? ctx5.

[4] 59w In channel pairs ctx5 ? ctx6 and ctx6 ? ctx5 the coupling level quickly restored after the initial decrease to values, which are significantly higher than the baseline level. These pairs were marked as ''primary". In channel pairs ctx4 ? ctx6 and ctx6 ? ctx4 the coupling restored only in the last 3 s. These pairs were marked ''secondary".

[5] 101w As can be seen in Fig. 3 and Table 2, the earliest significant increases in PI 0 (t) between cortico-thalamic channel pairs started to occur in the interval from 2 to 1 s before SWD onset and was characterized by gradual increase preceding SWD onset: ctx4 ? ATN, ctx4 ? VPM, ctx4 ? cRTN, ctx5 ? PO, ctx5 ? VPM, ctx5 ? cRTN, ctx6 ? PO, ctx6 ? cRTN. A bit later, 1-0.5 s preceding SWD onset other pairs became involved except ctx4 and ctx5 ? rRTN as well as ctx6 ? VPM, in which no increase preceding SWD was found.

[6] 111w A significant decrease in PI 0 (t) following a significant preictal increase was found for the following corticothalamic pairs: ctx4 ? PO, ctx4 ? VPM, ctx4 ? ATN, ctx5 ? PO, ctx5 ? VPM, ctx5 ? ATN, ctx6 ? PO, ctx6 ? rRTN, ctx6 ? ATN. A significant ictal decrease was also found for channel pair ctx4 ? rRTN, which did not show a preictal increase in PI 0 (t). The decrease always occurred in the first second after SWD onset and this decrease did not last longer than 1.5 s. There was also a combination of pairs for which the preictal increase was not followed by decoupling: cortex ? cRTN.

[7] 237w Instead, these cortex ? cRTN pairs showed an immediate (at SWD onset) or less than 1 s after SWD onset a primary increase in coupling which persisted until the end of the SWDs. Similarly, also cortex ? ATN channel pairs regained an increased coupling, following the initial drop at SWD onset, and the elevated increase persisted until SWD offset, i.e. it is proposed that for SWD maintenance the cortex starts to influence the Fig. 2. Dynamics of adapted nonlinear Granger causalities for intracortical channel pairs. Y-axis: PI, normalized to baseline level (10-3 s before onset). time, the moment of SWD onset is considered to be at t = 0. Prediction improvement (PI) was averaged per rat (16 rats). Black vertical lines indicate the seizure onset and offset, gray vertical lines indicate the length of moving window, in which Granger analysis was performed. Black points indicate values significantly larger than zero (baseline PI from 10 to 3 s prior to SWD is treated as a zero level) based on Student t-test, red points indicate values significantly lower than the preictal maximum, gray points -all others. In all cases only if there are at least three significant values in a row, they are considered as truly significant in order to avoid false positive due to repetitive testing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

[8] 55w ATN and cRTN as early as possible. Secondary increases were established for ctx5 ? PO and cortex ? rRTN; i.e. for these pairs an additional late occurring increase in coupling strength (varying between 3.5 till 0.5 s before the end of SWD) was noted following their preictal increase and drop in activity at SWD onset.

[9] 220w The dynamics regarding thalamo-cortical coupling are presented in Fig. 4 and Table 2. The earliest preictal increases in PI 0 (t) started to occur in the interval 2-1 s before SWD onset: ATN ? ctx5, ATN ? ctx6, PO ? ctx4, PO ? ctx5, PO ? ctx6, VPM ? ctx4, VPM ? ctx5. Fig. 3. Dynamics of adapted nonlinear Granger causalities for cortico-thalamic channel pairs. Y-axis: PI, normalized to baseline level (10 to 3 s before onset). X-axis: time, the moment of SWD onset is considered to be at t = 0. Prediction improvement (PI) was averaged per rat (16 rats). Black vertical lines indicate the seizure onset and offset, gray vertical lines indicate the length of moving window, in which Granger analysis was performed. Black points indicate values significantly larger than zero (baseline PI from 10 to 3 s prior to SWD is treated as a zero level) based on Student t-test, red points indicate values significantly lower than the preictal maximum, gray points -all others. In all cases only if there are at least three significant values in a row, they are considered as truly significant in order to avoid false positive due to repetitive testing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

[10] 67w Next, all other interactions from thalamus to cortex increased except rRTN ? cortex. Interestingly, many of these channel pairs which showed a significant increase were the same ones which also showed an increase in cortico-thalamic coupling. This was also the case for the interaction between cortex and cRTN: it was first are unidirectional (cortex ? cRTN) and it became bidirectional only in 1 s before SWD onset.

[11] 129w An immediate (at SWD onset) and significant decrease in PI 0 (t) lasting up to 1.75 s after the SWD onset, following a significant preictal increase which was found for: ATN ? ctx4, ATN ? ctx5, cRTN ? ctx4, cRTN ? ctx5, PO ? ctx4, PO ? ctx5, VPM ? cortex. A significant ictal decrease was also found for some channel pairs, which did not show a preictal increase in PI 0 (t): rRTN ? ctx4, rRTN ? ctx5. A decrease in coupling was found from the thalamus to the 4th and 5th somatosensory cortex layers, rather than to the 6th layer. Layers 4 and 5 are receiving inputs from the thalamus, while layer 6 drives the VPM and layer 5 drives higher order nuclei, among others the PO.

[12] 28w In thalamo-cortical channel pairs there was no primary increase. Secondary increases were found for ATN ? cortex and cRTN ? cortex. Notice that also these interactions became bidirectional.

[13] 51w Intrathalamic coupling changes (the data are presented in Fig. 5 and Table 2) were also noticed 2-1 s before the SWD onset: they were represented by bidirectional increase in two channel pairs (cRTN M Po and VPM M Po) and a unidirectional increase between ATN ? VPM, and rRTN ? PO.

[14] 63w The significant and immediate decrease in PI 0 (t) following a significant preictal increase was found for the following intrathalamic pairs VPM M PO, ATN ? VPM, and rRTN ? PO. A significant ictal decrease was also found for some channel pairs, which did not show a preictal increase in PI 0 (t): ATN ? PO, PO ? rRTN, and rRTN ? VPM.

[15] 39w Primary increases can be seen for following channel pairs: ATN M cRTN and cRTN M PO, somewhat later for cRTN M VPM and cRTN ? rRTN, suggesting that bidirectional intrathalamic interactions are important in the initial SWD maintenance phase.

[16] 58w Secondary increases in intrathalamic network were found for ATN M rRTN (first unidirectional, later bidirectional), VPM ? ATN, and rRTN ? cRTN. This implies that new interactions (bidirectional from the start) emerged which lasted till the end of the SWD, and, that in addition, the interaction between rRTN and cRTN became bidirectional toward the end of the SWD.

[17] 61w The method considers two signals, X and Y that were recorded simultaneously from two brain areas and that were used for constructing two time series, fx n g N n¼1 and fy n g N n¼1 ,correspondingly at the n-th time point. The study of causal (driving) interactions between X and Y with the aid of Granger causality includes three steps.

[18] 27w First, a univariate predictive model was constructed based on the one-channel raw data fx n g N n¼1 in the form of model map (Eq. ( 1)):

[19] 264w where x 0 nþs is the predicted value corresponding to the measured value x n+s (Fig. 4); f is approximating function as polynomials with order P; xn ¼ ð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), which is a classical approach to transpose time series in phase space, i.e. to obtain the high-dimensional state vector fx n g N0 n¼1 from the scalar time series fx n g N n¼1 for each time point, where N 0 ¼ N À s À maxððD s À 1Þl; l T Þ is an efficient time series length; s 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; D s is the embedding dimension that is actually the number of components in a state vector (Kugiumtzis, 1996); l is time delay (or lag), i.e. time interval between EEG values is used to construct the state vector; l T is the additional lag that takes into account the value of the experimental data delayed from the indicial (predicted) time point with a period of T (Sysoeva and Sysoev, 2012). The difference between the predicted x 0 nþs and observed x nþs values is the prediction error (in Fig. 4(b) vertical distance between x nþs and x 0 nþs ). Model coefficients were selected using least square estimates (Legendre, 1805), i.e. by minimizing the squared prediction error, e 2 s (Eq. ( 2)):

[20] 56w where r 2 is the empirical dispersion of time series fx n g N n¼1 , n s ¼ max ðD s À 1Þl; l T ð Þ þ 1. Second, the bivariate model (Eq. ( 2)) was constructed based on both time series fx n g N n¼1 and fy n g N n¼1 :

[21] 24w x 00 nþs ¼ g x n ; x nÀl ; . . . x nÀðDsÀ1Þl ; y n ; . . . ;

[22] 80w x 00 nþs is the predicted value corresponding to the measured value x nþs , received two time series fx n g N n¼1 and fy n g N n¼1 ; D a is dimension of the state vector ỹn ¼ ðy n ; y nÀl ; . . . y nÀðDaÀ1Þl Þ reconstructed from the scalar time series fy n g N n¼1 in (Eq. ( 3)), so the total dimension of the bivariate model can be computed as

[23] 258w Y is stated to drive X if the value x 00 nþs predicted with the bivariate model is closer to the measured one x nþs than the value predicted with the univariate model x 0 nþs (in average). This principle is illustrated in Fig. 7, where the vertical distance between x 00 nþs (circle) and x nþs (triangle) is smaller than the vertical distance between x 0 nþs (asterisk) and x n+s (triangle): in such a case at the time point n the bivariate prediction error e 2 j is smaller than the univariate one e 2 s . The prediction length value in the model (Eqs. ( 1) and ( 2)) was set to s = T/4, where T is duration of one characteristic period in a signal (Sysoeva et al., 2014). The chosen value of s provided the best compromise between sensitivity of the presented method (capability to detect the actual coupling) and specificity (minimum of false positive decisions about coupling). Bayesian information criterion (BIC) was used in order to determine the optimal values of P and D s (Schwarz, 1978). The value of lag l was optimized based on the fact that, in systems with characteristic period T, the best specificity could be obtained with the lag varying between T/12 and T/3 (Kornilov et al., 2014), and the values l = 1 and l = n _ s T/2 (where n is an arbitrary natural number) must be avoided, because of a very high probability of false positive results. Therefore we set l = T/10.

[24] 29w Third, the value of prediction improvement PI was computed with Eq. ( 4), and it is considered as the most important measurable characteristic of the adapted Granger causality method.

[25] 92w The situation when e 2 j ¼ e 2 s suggests that the data from the second EEG channel Y, i.e. time series fy n g N n¼1 , do not improve a prediction of fx n g N n¼1 . In other words, Y does not drive X. The situation when e 2 s > 0 and e 2 j ! 0, so PI ? 1 suggests that the data from the second time series fy n g N n¼1 significantly improve prediction of the first one, so Y drives X.

[26] 9w (Accepted 18 November 2015) (Available online 2 December 2015)

[1] 417w Processes involved in the occurrence, maintenance, and termination of SWD might be invisible on the raw EEG records, since they might not take place in a particular brain structure but between different structures. However, they can be detected using time-resolved network analysis methods such as coupling analysis applied to multichannel time series, such as time-variant adapted nonlinear GC The nonlinear adapted GC method has a high sensitivity (a capability to detect the actual coupling) and a high specificity (the number of false positive results is small) (Sysoeva et al., 2014). This high sensitivity is achieved by means of the nonlinearity of Fig. 4. Dynamics of adapted nonlinear Granger causalities for thalamo-cortical channel pairs. Y-axis: PI, normalized to baseline level (10 to 3 s before onset). X-axis: time, the moment of SWD onset is considered to be at t = 0. Prediction improvement (PI) was averaged per rat (16 rats). Black vertical lines indicate the seizure onset and offset, gray vertical lines indicate the length of moving window, in which Granger analysis was performed. Black points indicate values significantly larger than zero (baseline PI from 10 to 3 s prior to SWD is treated as a zero level) based on Student t-test, red points indicate values significantly lower than the preictal maximum, gray points -all others. In all cases only if there are at least three significant values in a row, they are considered as truly significant in order to avoid false positive due to repetitive testing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) the method. It has been shown in a number of works that the use of nonlinear approximating functions makes it possible to detect couplings, which cannot be detected using linear ones (Le van Quyen et al., 1999;Lehnertz, 1999;Lopes da Silva et al., 2003;Chen et al., 2004;Sysoev et al., 2010). In fact, there are also other reasons for the application of nonlinear methods in the analyses of LFP: networks that consist of a combination of excitatory and inhibitory cells show nonlinear behavior (Buzsa´ki, 2006) and the cortex and thalamus consists of combinations of excitatory and inhibitory cells. The application of a specially selected prediction length makes it possible to improve the specificity of the method (Sysoeva et al., 2014). The use of nonuniform embedding allows a reduction in the number of model coefficients and improves the temporal resolution of the method, using a short time window (only four main oscillations).

[2] 152w Bidirectional, unidirectional and a lack of coupling can be distinguished in most cases. In particular, unidirectional coupling can be detected even in the state that it is very close to phase synchronisation (for values of the phase synchronisation index (Allefeld and Kurths, 2004) up to 0.95 (Kornilov et al., 2014;Sysoev and Sysoeva, 2015)). The novel method is stable for synchrony and synchronisation, signal/noise ratio (amplitude), mediated (indirect) couplings (Sysoev and Sysoeva, 2015). Complete synchronization, however, prevents the detection of the coupling directionality, since states of interacting subsystems become undistinguishable. And that is often seen in many channel pairs during SWD: many of them show low and nonsignificant values during SWD when the cortico-thalamic network is in a highly synchronized state (McCormick and Contreras, 2001;Meeren et al., 2002). Another restriction of the method is that it cannot describe fast transient processes, this phenomenon of ''Ears" has been described elsewhere (Sysoev and Sysoeva, 2015).

[3] 252w The application of a new nonlinear time-variant adapted GC method describing intracortical, corticothalamic, thalamo-cortical and intrathalamic communications before, during and after the occurrence of SWDs confirmed the early involvement of the cortex in the initiation of the SWDs, next it showed a completely new process: a drop in coupling in the first second after SWD onset between the cortical layers, and in a subset of cortico-thalamic and intrathalamic Fig. 5. Dynamics of adapted nonlinear Granger causalities for intrathalamic channel pairs. Y-axis: PI, normalized to baseline level (10 to 3 s before onset). X-axis: time, the moment of SWD onset is considered to be at t = 0. Prediction improvement (PI) was averaged per rat (16 rats). Black vertical lines indicate the seizure onset and offset, gray vertical lines indicate the length of moving window, in which Granger analysis was performed. Black points indicate values significantly larger than zero (baseline PI from 10 to 3 s prior to SWD is treated as a zero level) based on Student t-test, red points indicate values significantly lower than the preictal maximum, gray points -all others. In all cases only if there are at least three significant values in a row, they are considered as truly significant in order to avoid false positive due to repetitive testing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) channel pairs. It showed also a maintenance process that developed into the spontaneous abortion of the SWD.

[4] 77w The SWD initiation process was characterized by an early preictal (more than 2 s prior SWD onset) increase in coupling within the intracortical layers of the somatosensory cortex. This well agrees and expands our early findings achieved with frequency resolved GC (Lu¨ttjohann et al., 2014) and with the now widely accepted view that SWDs have a cortical origin and are preceded by precursor and network activity (Meeren et al., 2002;Polack et al., 2007;van Luijtelaar et al., 2011).

[5] 67w Then the thalamus became more involved in the network activity. The cortex drives all thalamic nuclei except the rRTN, while VPM, PO and ATN respond. Just before SWD onset, cRTN also begins to drive the cortex. So most cortico-thalamic and thalamo-cortical pairs begin to interact, but rRTN and cortex remain uncoupled. Intrathalamic pairs are faintly involved in the initiation: only 1/3 of them participate in this process.

[6] 199w The increase in coupling in almost all cortico-cortical and cortico-thalamo-cortical pairs directly before SWD onset well agrees with EEG data from another EEG dataset of the same absence model (Sysoeva et al., 2014). On the other hand, the newly developed nonlinear method seems to be more sensitive than the frequencyresolved linear version of GC (Lu¨ttjohann et al., 2014) since more channel pairs showed significant results while the same data set was used. The adapted nonlinear GC method gave also more positive results in comparison to the nonlinear association analysis (Lu¨ttjohann and van Luijtelaar, 2012). Different statistical approaches (repeated measures analyses of variance and one-sided paired t-post hoc tests versus 2-sided t-tests) and different control periods (a 540-ms duration segment of passive wakefulness 5 min remote from the SWD occurrence, versus a 7 s interval from 10 to 3 s before SWD onset) might be the reason for the higher sensitivity. Also, earlier a directional coupling was only inferred if PI 0 (t) values for the direction from X to Y were significantly higher than in the opposite direction; in the current study the coupling detection in directions from X to Y and from Y to X was performed independently.

[7] 43w In all, the data suggest that SWDs emerge from the cortex and that preictally, the cortico-thalamo-cortical network gradually expands and later involves also intrathalamic coupling, and that the gradual increasing coupling becomes more and more bidirectional until the moment at which SWDs emerge.

[8] 109w As can be seen in Figs. 2345, many channel pairs lost their directional coupling at SWD onset: the vast majority of the channel pairs with preictal increases in PI 0 (t) showed a temporary drop immediately after the SWD onset. PI 0 (t) returned to baseline values within 1-1.5 s. This decrease occurred in 57% (32/56) of the channel pair combinations; it was detected in the majority of the intracortical (4/6), cortico-thalamic (10/15) and thalamo-cortical (11/15) channel pairs, but only for a minority (7/20) of the intrathalamic pairs. Interestingly, a significant drop in PI 0 (t) was never noticed for channel pairs which directed their influence to the cRTN.

[9] 293w This decrease in PI 0 (t) immediately following the SWD onset suggests temporarily decoupling. Interestingly, this early decoupling was neither noticed with a linear frequency resolved GC method, nor with nonlinear association analyses (Lu¨ttjohann and van Luijtelaar, 2012;Lu¨ttjohann et al., 2014). Apparently, the lack of detecting with linear methods can be explained by a large synchrony between channel pairs, initiated preictally. Linear methods are often unable to distinguish between synchrony and synchronisation (Sysoev and Sysoeva, 2015). In order to establish whether this drop is genuine and not an artifact of the adapted GC, the mutual information function (MI) was calculated (Kraskov et al., 2004) for all channel pairs of all 16 rats of the original data set. The data are presented in Fig. 6. A decrease following SWD onset can be appreciated for the mutual information in a number of channel pairs, which also showed a drop in adapted nonlinear GC. This decrease is significant for six pairs: ctx4-ctx6, ctx4-ATN, ctx4-VPM, ctx4-cRTN, ctx5-VPM, VPM-cRTN. It has to be mentioned that the straightforward comparison of GC and MI is less meaningful, since MI is neither able to show direction of coupling nor its causality. It only demonstrates that the amount of common information in the signals dropped during the first second of a SWD and then restored. Actually, the GC can reveal the changes in coupling (the reason), while the MI only can show, how these changes reflect network dynamics (the result). This might also explain why the drop in MI at SWD onset occurs a bit later than in GC (first cause, followed by the effect). Assuming that this is indeed true, that it can be concluded that the outcomes of GC and MI analyses regarding the presence of decoupling confirm each other.

[10] 141w This decoupling is not easy to explain in neurophysiologic terms without concomitant simultaneous intra or extra cellular recordings in different regions. This decrease of PI 0 (t) was most often, although not exclusively, present after a gradual increased preictal coupling. It is thinkable that the preictal coupling increase pushes the network over a threshold and next synchronous oscillations in the form of SWD emerge. After SWD are initiated, and assuming that synchrony of oscillations is metabolically cheap (Buzsa´ki, 2006), it is proposed that only a few elements of the network need to interact (drive) in order to maintain oscillating for a while. The only remaining driving force in the first second of SWD is the drive from the deep cortical layers to the cRTN. This might be interpreted as that this cortical drive to the cRTN is sufficient for SWD maintaining.

[11] 97w From the point of view of oscillation theory, the initiation of SWD can be considered as a gradual increasing push, when a lot of energy is injected in the network. This energy was enough to keep the oscillatory activity for the 1st second. Since this energy might dissipate due to insufficient interactions, the oscillation frequency of SWDs decreases during the first second (see time frequency plots at Fig. 1 for most of the channels); frequency modulation of SWD in patients and in WAG/Rij has been reported previously (Bosnyakova et al., 2007). Only then another process -maintenance, starts.

[12] 102w The SWD maintenance process starts from cortical (ctx5, ctx6) influence to cRTN less than 1 s after the onset. Here, cRTN works as a transmitter, immediately starting to drive PO and ATN (also <1 s after onset). A little bit later cRTN becomes coupled with thalamic nuclei: PO, VPM, ATN (all bidirectional). Also intracortical communications appear and all cortical layers begin to drive ATN. Note, that there is no driving from thalamus to cortex during the first 3 s of seizure. In general, an early role of the cortex and a more passive involvement of thalamus in SWD activity maintenance were established.

[13] 28w Before the SWD offset the cortex is bidirectional coupled with cRTN and ATN, also cRTN is coupled bidirectional with all thalamic nuclei, playing central role in SWD maintenance.

[14] 76w The clear involvement of the ATN in the maintenance process of SWD was not revealed previously. The ATN, which is reciprocally connected with the rRTN and not directly with the cortex, seems only passively involved considering that lesions of the ATN in GAERS did not affect the incidence of SWD (Marescaux et al., 1992) and a passive role of the ATN is SWD occurrence has been proposed using nonlinear association analysis (Lu¨ttjohann and van Luijtelaar, 2012).

[15] 345w It is proposed that SWD termination can occur either due to the end of maintenance process, or due to that there is a third special process, being manifested by a relatively short increase in coupling just before seizure offset. In the first case a permanent increase in coupling toward the end of a SWD is present, perhaps necessary to maintain the SWD, which can be seen in the following pairs: ctx5 M ctx6; cortex ? cRTN; cortex ? ATN; cRTN M ATN; cRTN M PO; cRTN M rRTN; cRTN M VPM. To maintain such a high level of coupling is physically very energy consuming, i.e. a lot of energy is necessary. Therefore this increase cannot last for a long time and neurons become less hyperpolarized (Polack et al., 2007), interspike frequency may diminish toward the end of SWD (Bosnyakova et al., 2007), the amplitude of the spikes of the SWDs diminishes, and the spikes become less sharp, depression of coupling occurs, the various previously coupled pairs desynchronize and SWDs are terminated. Our findings suggest that the maintenance process needs intracortical interactions, Fig. 6. Function of mutual information calculated in moving window of 0.5 s length with 0.125-s shift between windows from 10 s before to 4.5 s after the seizure onset, and averaged over all seizures for 16 rats with all channels available. The gray, black and red points have the same meaning as for Fig. 2345. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) cortico-thalamic interactions (mainly driving from the cortex to the cRTN and ATN), and also bidirectional coupling of the cRTN with all thalamic nuclei. The second case can be seen in all rRTN M cortex pairs. During the second half of seizure the cortex drives rRTN. When the rRTN responds to this coupling, the seizure terminates, most likely via a negative feedback loop. Interestingly, this contribution of these channel pairs in SWD termination is unique since rRTN driving the cortex did not occur in any other stage.

[16] 97w van Luijtelaar G, Sitnikova E, Lu¨ttjohann A (2011) On the origin and suddenness of absences in genetic absence models. Clin EEG Neurosci 42(2):83-97. Westmijse I, Ossenblok P, Gunning B, van Luijtelaar G (2009) Onset and propagation of spike and slow wave discharges in human absence epilepsy: a MEG study. Epilepsia 50(12):2538-2548. Vergnes M, Marescaux C, Depaulis A, Micheletti G, Warter J (1987) Spontaneous spike and wave discharges in thalamus and cortex in a rat model of genetic petit mal-like seizures. Exp Neurol 96:127-136. Wang CW (2007) Nonlinear phenomena research perspectives. New York: Nova Science Publishers. pp 7-53.

[1] 109w Sixteen male WAG/Rij rats, 6-9 months of age were used as experimental subjects. They were born and raised at the department of Biological Psychology, Donders Centre for Cognition, Radboud University Nijmegen, The Netherlands. Prior to surgery rats were housed in pairs (High MakrolonÒ cages with Enviro DriÒ bedding material and cage enrichment) with free access to food and water and were kept at a 12-12 h light-dark cycle (light off at 8.30 AM). After surgery rats were housed individually. The experiment was approved by the Ethics Committee on Animal Experimentation of Radboud University Nijmegen (RU-DEC). Efforts were made to keep the discomfort for the animals as minimal as possible.

[2] 234w Implantation of the LFP recording electrodes was done in a stereotactic frame under isoflurane anesthesia. At the start of surgery, rats received a subcutaneous injection of the analgesic RimadylÒ and an intramuscular injection of atropine to prevent excessive salivary production. Body temperature was controlled and conserved via a heating pad. The local anesthetic Lidocaine was used on the incision points. Holes were drilled into the skull on top of the right hemisphere for the insertion of recording electrodes at the following positions: Somatosensory cortex: A/P = 0.0, M/L = À4.6 depth = À2.8 (layer 4), À3.1 (layer 5), À 3.6 (layer 6); anterior thalamus: A/ P = À1.4 M/L = À1, depth = À6.2; rRTN: A/P = À1.4, M/L = À1.9, depth = À6.6; PO: A/P = À3.6, M/L = À2, depth = À5.4; VPM: A/P = À4.16, M/L = À2.8, depth = À6 and cRTN: A/P = À3.1, M/L = À3.5, depth = À6.6. All coordinates were determined relative to bregma according to the rat-brain atlas of Paxinos and Watson, 2006. Electrode wires, assembled in a self-constructed electrode system (Lu¨ttjohann and van Luijtelaar, 2012) were simultaneously inserted into the brain. Ground and reference electrodes were positioned epidurally on top of the cerebellum. The electrode assembly was fixed to the skull via dental cement. Postoperative analgesic RimadylÒ (24 and 48 h after surgery) was administered, and rats were allowed to recover for two weeks.

[3] 130w Only data from brain structures with a histologically verified proper electrode position were included in the statistical analysis (Table 1). To facilitate the finding of the location of the tip of the recording electrodes, a direct current (9 V, 25 lA, 10 s duration) was passed through each electrode in the deeply anesthetized rat at the end of the experiment. Next rats were perfused with a potassiumferrocyanide-formaldehyde-phosphate solution, coloring these lesions at the end of each electrode tip. Brains were fixed in a 30% sucrose solution, 0.1-ml PBS, cut in 40-lm coronal slices with the aid of a microtome, and stained with Cresyl Violet. Only electrodes for which the midpoint of the small lesion was located within the target structure were considered as properly implanted and included in statistical analysis.

[1] 137w Two weeks after surgery rats were placed individually in a 20 Â 35 Â 25 inch Plexiglas registration boxes and connected to the recording leads for multichannel LFP recordings. These were attached to a swivel-contact, which allowed recording in freely moving animals. The LFP signals were amplified with a physiological amplifier (TD 90087, Radboud University Nijmegen, Electronic Research Group), filtered by a band pass filter with cutoff points at 1 (high pass) and 100 (low pass) and a 50 Hz Notch filter, and digitized with a constant sample rate of 2048 Hz on a WINDAQ recording system (DATAQ-Instruments). The movements of rats were registered by means of a Passive Infrared Registration system (PIR, RK2000DPC LuNAR PR Ceiling Mount, Rokonet). Each rat was recorded for a period of 4 h during the dark phase of the light-dark cycle.

[2] 46w In order to estimate time-dependent changes of coupling characteristics between different brain areas and especially with the onset and offset of SWD, a specially designed adapted GC approach (Sysoeva et al., 2014) was used based on principles of a time-variant GC method (Hesse et al., 2003).

[3] 213w The main advantage of GC is that it is less dependent on the amount of experimental data as most other methods. This is due to parameterization: one has to estimate a small number of coefficients of univariate and bivariate models, rather than multidimensional distributions as for transfer entropy. This gives the opportunity to analyze nonstationary data in a moving window. However, wrong parameterization dramatically suppresses the efficiency of the method, giving either a lot of false positive results, or failures to find actually present coupling. The main problems could be: insufficient sampling rate (Smirnov and Bezruchko, 2012), inadequate consideration of signal spectral characteristics when the model parameters for prediction length and time lag have to be chosen (Sysoeva and Sysoev, 2012;Kornilov et al., 2014), insufficient number of nonlinear terms (Sysoev et al., 2010) or inadequate chosen nonlinear functions (Chen et al., 2004), all leading to wrong model dimensions. To solve these problems, a special structure of an empirical model for absence seizures was developed (Sysoeva and Sysoev, 2012) and this was applied in a two (frontal cortex and VPM) channel EEG recording for the analysis of coupling in WAG/Rij rats (Sysoeva et al., 2014). Now the same empirical model was applied here on the data set published in (Lu¨ttjohann and van Luijtelaar, 2012).

[4] 239w Calculations of prediction improvement (PI) were performed in a moving window of 500 ms (1024 data points) that was shifted in time by 125 ms (256 data points), since SWD changes in connectivity are known to occur within timeframes as short as 0.5 s (Meeren et al., 2002). The method is described in more detail in Appendix A. Please note, that PI can only be used for detecting coupling presence or changes in coupling. Its absolute value cannot be interpreted equivocally. Also analysis in two directions, e.g. VPM ? ctx4 and ctx4 ? VPM, is completely independent in all cases. This means that one does not have to compare the PI values in opposite directions to understand the coupling directionality. e.g. let PI VPM!ctx4 ¼ 0:4, if this value is significantly different (e.g. with p-value <0.05) from baseline score as determined for a channel pair in a specific direction, we accept that there is coupling in this specific direction, in our example from VPM to ctx4. Also let PI ctx4!VPM ¼ 0:5, than the base-line level of PI (the other one) is also calculated in the other direction, so ctx4 ? VPM. If PI ctx4!VPM is significantly different from baseline in the direction ctx4 ? VPM, the coupling has become significant as well, so the coupling is bidirectional. However, if PI ctx4!VPM differs from its baseline level insignificantly, that the coupling is unidirectional (only VPM ? ctx4), even though

[5] 197w Here, 10 fragments containing an SWD were selected from each animal (155 fragments in total, in one rat only five SWDs were included because of technical reasons). This number of fragments was chosen since the intraindividual variation in morphology in SWD is small compared to its interindividual variation and this number was the minimum of SWDs with an acceptable length of at least 5 s. Each considered LFP fragment containing SWD (see example at Fig. 1) was split into two parts: 2. Summary of interactions in cortico-thalamo-cortical network for all phases of seizure development. Channels: ctx4-6 are somatosensory cortex layers, ATN -anterior thalamic nucleus, PO -posterior thalamic nucleus, VPM -ventro-posteromedial thalamic nucleus, rRTN and cRTN are the rostral and caudal parts of reticular thalamic nucleus. For the initiation process all channel pairs are divided into four groups: (1) >2 s before SWD onset, (2) in the interval 2-1 s before SWD onset, (3) less than 1 s before SWD onset, (4) no initiation. For the maintenance process two types are shown: (1) primary, when the process starts just or little bit after the decoupling and (2) and secondary, when it starts at the second half of seizure

[6] 567w Channel pair Initiation Decoupling Maintenance Termination Cortico-cortical channel pairs ctx4 ? ctx5 >2 s before SWD onset + ctx4 ? ctx6 >2 s before SWD onset Secondary, >2 s before SWD offset ctx5 ? ctx4 61 s before SWD onset + ctx5 ? ctx6 >2 s before SWD onset Primary, 1-2 s after SWD onset ctx6 ? ctx4 >2 s before SWD onset + Secondary, >2 s before SWD offset ctx6 ? ctx5 >2 s before SWD onset + Primary, 1-2 s after SWD onset Cortico-thalamic channel pairs ctx4 ? ATN 2-1 s before SWD onset + Primary, >2 s after SWD onset ctx4 ? cRTN 2-1 s before SWD onset Primary, 1-2 s after SWD onset ctx4 ? PO 61 s before SWD onset + ctx4 ? rRTN + Secondary, >2 s before SWD offset ctx4 ? VPM 2-1 s before SWD onset + ctx5 ? ATN 61 s before SWD onset + Primary, 1-2 s after SWD onset ctx5 ? cRTN 2-1 s before SWD onset Primary, 61 s after SWD onset ctx5 ? PO 2-1 s before SWD onset + Secondary, 2-1 s before SWD offset ctx5 ? rRTN Secondary, >2 s before SWD offset ctx5 ? VPM 2-1 s before SWD onset + ctx6 ? ATN 61 s before SWD onset + Primary, 1-2 s after SWD onset ctx6 ? cRTN 2-1 s before SWD onset Primary, 61 s after SWD onset ctx6 ? PO 2-1 s before SWD onset + ctx6 ? rRTN 61 s before SWD onset + Secondary, >2 s before SWD offset ctx6 ? VPM Thalamic-cortical channel pairs ATN ? ctx4 61 s before SWD onset + Secondary, >2 s before SWD offset ATN ? ctx5 2-1 s before SWD onset + Secondary, 2-1 s before SWD offset ATN ? ctx6 2-1 s before SWD onset Secondary, >2 s before SWD offset cRTN ? ctx4 61 s before SWD onset + Secondary, >2 s before SWD offset cRTN ? ctx5 61 s before SWD onset + Secondary, 2-1 s before SWD offset cRTN ? ctx6 61 s before SWD onset Secondary, >2 s before SWD offset PO ? ctx4 2-1 s before SWD onset + PO ? ctx5 2-1 s before SWD onset + PO ? ctx6 2-1 s before SWD onset rRTN ? ctx4 + + rRTN ? ctx5 + + rRTN ? ctx6 Secondary, >1 s before SWD offset + VPM ? ctx4 2-1 s before SWD onset + VPM ? ctx5 2-1 s before SWD onset + VPM ? ctx6 61 s before SWD onset + Intrathalamic channel pairs ATN ? cRTN Primary, 1-2 s after SWD onset ATN ? PO + ATN ? rRTN Secondary, 2-1 s before SWD offset ATN ? VPM 2-1 s before SWD onset + cRTN ? ATN Primary, 61 s after SWD onset cRTN ? PO 1-0.5 s before SWD onset Primary, 61 s after SWD onset cRTN ? rRTN Primary, >2 s after SWD onset cRTN ? VPM Primary, 1-2 s after SWD onset PO ? ATN PO ? cRTN 61 s before SWD onset Primary, 1-2 s after SWD onset PO ? rRTN + PO ? VPM 2-1 s before SWD onset + rRTN ? ATN Secondary, >2 s before SWD offset rRTN ? cRTN Secondary, >2 s before SWD offset rRTN ? PO 2-1 s before SWD onset + rRTN ? VPM + (continued on next page)

[7] 41w 1. preictal (10 s prior to SWD onset) to ictal transition (3 s after SWD onset), 2. ictal (3 s prior to SWD offset) to postictal transition (10 s after SWD offset) with possible small overlapping for 17 of 160 segments.

[8] 276w SWD onset was determined upon visual inspection of the LFP recordings by a trained electrophysiologist. As reported previously (Lu¨ttjohann and van Luijtelaar, 2012;Lu¨ttjohann et al., 2013Lu¨ttjohann et al., , 2014) ) the onset of a SWD was defined by the presence of a sharp spike, with an amplitude of at least two times of the background, which is present in both cortical and thalamic channels and is followed by rhythmic SWD activity. Likewise the offset of SWD was represented by the last, high-amplitude spike, which was present in all cortical and thalamic channels. In the most of the cases (103 of 155), the rats were completely immobile in all 3 s before SWD onset and during the SWDs. In most other cases (41 of 155) they were immobile with small head movement. Only in 11 cases of 155 (7%) rats were actively moving preceding the SWD. This was controlled by inspecting the PIR channel (an infrared movement detector, which gave an analog signal in parallel with the LFP recordings), indicating movement activity. Resulting dependencies of PI on time were averaged across all seizures in each animal, matching start and ending moments of seizures. Then for each averaged dependency PI(t) the background level PI bg was established as an average PI over 7 s time interval (baseline period, from 10 to 3 s before SWD start). This period was completely devoid of SWDs. Normalized dependencies were calculated as PI 0 ðtÞ ¼ PIðtÞ À PI bg . The value of PI 0 (t) = 0 corresponds to the baseline level, positive values of PI 0 correspond to larger coupling than in baseline and negative -to lower one.

[9] 136w Series of PI 0 (t) for each rat were averaged forming a sample of PI 0 values (from 10 to 16 for different channel pairs, depending on the availability in the LPF data set) for each time point. These data points were analyzed with Student t-tests to establish differences from zero with p-value <0.05. In order to correct for an increased chance of false positives and to reduce the chance of getting type I errors, only clusters of minimum three subsequent significant timepoints were considered as representing as showing genuine changes in coupling. Therefore, a SWD transition (initiation, maintenance, abortion) process is only considered if there are not less than three significantly subsequent points in PI 0 (t) different from zero (i.e. from baseline) with an exception of the areas between the gray and black lines.

[10] 104w It was noticed that the initiation process was followed by a decrease in PI 0 (t) (after the gray line) in many channel pairs. This drop in PI 0 (t) was statistically evaluated with paired t-tests by taking the difference between the preictal PI maximum in the time window from 2 to 0.5 s prior SWD onset, the exact time point of this maximum may vary from one channel pair to another, and the PI minimum in the first 1.5 of SWD after SWD onset. Again, only a cluster of three subsequent significant timepoints in PI 0 (t) were considered as genuine significantly different.

[11] 150w Vertical black lines on PI 0 (t) plots indicate the seizure onset or offset timepoints, gray vertical lines indicate the length of the moving window, in which Granger analysis was performed. The results between the black and gray lines have to be considered as unsafe due to effect of a transition state (when the model that was constructed is partly based on the previous regime, partly on the next one) and will not be considered. This phenomenon has been studied on etalon oscillators and was named ''Ears" (Sysoev and Sysoeva, 2015). ''Ears" are an increase of PI, when the moving window covers the fast transition from the preictal to the ictal and from the ictal to the postictal phase. Therefore, significant points between the black and gray lines are not taken into account. Therefore, only timepoints from 1 to 0.5 s before SWD onset and before SWD offset were interpreted.