## Abstract

Inferring strength and direction of interactions from electroencephalographic (EEG) recordings is of crucial importance to improve our understanding of dynamical interdependencies underlying various physiological and pathophysiological conditions in the human epileptic brain. We here use approaches from symbolic analysis to investigate—in a time-resolved manner—weighted and directed, short- to long-ranged interactions between various brain regions constituting the epileptic network. Our observations point to complex spatial–temporal interdependencies underlying the epileptic process and their role in the generation of epileptic seizures, despite the massive reduction of the complex information content of multi-day, multi-channel EEG recordings through symbolization. We discuss limitations and potential future improvements of this approach.

## 1. Introduction

Epilepsy is one of the most common neurological disorders, second only to stroke, that affects approximately 1% of the world's population [1,2]. In about 25% of patients, epileptic seizures—the cardinal symptom of epilepsy—cannot be controlled by any available therapies [3–7]. In order to enhance the quality of life for epilepsy patients, there is thus a great need for improved therapeutic possibilities. An epileptic seizure is defined as ‘a transient occurrence of signs and/or symptoms due to abnormal, excessive or synchronous neuronal activity in the brain’ [8,9]. Epileptic seizures can be divided into two main categories: focal seizures, which appear to originate from a circumscribed brain region (the so-called seizure-onset zone (SOZ) [10]) and which may or may not remain restricted to this region, and generalized onset seizures, which appear to engage almost the entire brain. These concepts of focal and generalized seizures, however, are being challenged by increasing evidence of seizure onset within a network of brain regions (the so-called epileptic network) which led to a new approach to classification of seizures and epilepsies [11,12]. In an epileptic network, all its constituents can contribute to the generation, maintenance, spread and termination of seizures as well as to the many pathophysiological phenomena seen during the seizure-free interval [13–17]. The concept of an epileptic network may also explain the apparent contradictory findings of seizure precursors that are not confined to the SOZ or its immediate surroundings but can be observed in remote or even contralateral brain regions [18–27].

An improved understanding of the epileptic network and its complex dynamics can be achieved with time-series analysis techniques that aim at characterizing interactions and their properties, namely strength and direction [28–35]. Over the past years, a number of such bivariate analysis techniques have been proposed, ranging from linear to nonlinear ones. Comprehensive overviews concerning these techniques and their applications in diverse fields can be found in references [36–43]. Here, we concentrate on information theoretic approaches [44], and particularly on those that use symbol sequences [45–47] derived from empirical time series to characterize interactions. For the strength of interaction, we will consider the order parameter proposed by Liu [48] and for the direction of interaction, the symbolic transfer entropy [49,50]. Both approaches allow us for a robust and computationally fast quantification of the strength and the preferred direction of information flow between time series from observed data.

Previous studies [49,50] that investigated directed interactions in the human epileptic brain provided evidence that the aforementioned approaches allow one to reliably identify the hemisphere containing the SOZ without observing actual seizure activity. In these studies, however, analyses were restricted to data from patients with mesial temporal lobe epilepsies and to recordings from within the mesial temporal lobes. Here, we extend these studies in several aspects. First, we investigate multi-day (on average about 4 days), multi-channel electroencephalographic data recorded invasively from 11 patients that suffered from seizures with different anatomical onset locations. Second, we investigate—in a time-resolved manner—changes in interdependencies between all sampled brain regions, thus taking into account interactions for a variety of physiological and pathophysiological conditions. Third, we estimate *both* strength and direction of interactions in order to effectively distinguish the various coupling regimes (uncoupled, weak to strong couplings) that may be observable in the human epileptic brain and in order to avoid misinterpretations [31,50]. Fourth, we address the question of who is driving preferentially whom in large-scale epileptic brain networks and whether consistent changes in the directed interactions can be identified when approaching epileptic seizures.

## 2. Strength and direction of interaction from symbolic analysis

In the following, we denote with *x*_{i}:=*x*(*i*Δ*t*) and *y*_{i}:=*y*(*i*Δ*t*), *i*=1,…,*N*, time series of observables of systems *X* and *Y* recorded simultaneously with sampling interval Δ*t*. We make use of the seminal work of Bandt & Pompe [46] and derive symbol sequences from reordering the amplitude values of time series. Let *l* and *m* denote the embedding delay and embedding dimension, which were chosen appropriately for symbolization. Then, *m* amplitude values *s*_{i}={*x*_{i},*x*_{i+l},…,*x*_{i+l(m−1)}} for a given, but arbitrary timestep *i* are arranged in ascending order {*x*_{i+l}(*k*_{i1}−1)≤*x*_{i+l(ki2−1)}≤⋯≤*x*_{i+l(kim−1)}} with rank *k*_{ii}. In the case of equal amplitude values, the rearrangement is carried out according to the associated index *k*, i.e. for *x*(*i*+(*k*_{i1}−1)*l*)=*x*(*i*+(*k*_{i2}−1)*l*), we write *x*(*i*+(*k*_{i1}−1)*l*)≤*x*(*i*+(*k*_{i2}−1)*l*) if *k*_{i1}<*k*_{i2}. This ensures that every *s*_{i} is uniquely mapped onto one of the *m*! possible permutations, and a permutation symbol is defined as
2.1
The symbol sequence for system *Y* is defined in complete analogy.

A robust method to estimate the strength of interaction from symbol sequences was proposed by Liu [48]. It is based on consistent changing tendencies of temporal permutation entropies , where the sum runs over all possible symbols. *H*_{Y}(*w*_{η}) is defined in complete analogy. In order to estimate these changing tendencies, symbol sequences and , *i*=1,…,*N* are split into *N*_{η} segments *w*_{η} with *η*=1,…,*N*_{η}, consisting of *N*_{w} datapoints, and permutation entropies *H*_{X} and *H*_{Y} are estimated for each segment *w*_{η}. *H*_{X}(*w*_{η}) and *H*_{Y}(*w*_{η}) can be expected to be similar if some functional relationship—in the sense of generalized synchronization [51]—exists between states of *X* and *Y* . The changing tendency for a time series from system *X* can be quantified with
2.2
and in complete analogy for a time series *S*_{Y} from system *Y* . A measure for the strength of interaction between systems *X* and *Y* can then be defined as
2.3
This index will be around zero for independent time series and close to unity for generalized synchronization. Note that might attain slightly negative values.

Symbolic transfer entropy (STE) [49,50] has been proposed as a computational fast and robust method to quantify the direction of interaction between coupled systems. It is based on Schreiber's transfer entropy [52] that allows one to distinguish effectively driving and responding elements and to detect asymmetry in the interaction of systems. STE overcomes some of the limitations of previous techniques that aim at estimating transfer entropy from time series, such as requiring fine tuning of parameters and being highly sensitive to noise contributions [53–55]. More specifically, STE makes use of relative frequencies of symbols to estimate the transition probabilities required to quantify the deviation from the generalized Markov property, *p*(*x*_{i}|*x*_{i−1},*y*_{i−1})=*p*(*x*_{i}|*x*_{i−1}), where *p*(⋅|⋅) denotes the conditional transition probability density. With transfer entropy, this deviation is formulated as a Kullback–Leibler entropy between *p*(*x*_{i}|*x*_{i−1},*y*_{i−1}) and *p*(*x*_{i}|*x*_{i−1}).

With given symbol sequences and , symbolic transfer entropy is then defined as
2.4
where the sum runs over all symbols. is defined in complete analogy, and the directionality index allows one to quantify the preferred direction of information flow between systems *X* and *Y* . attains positive values for unidirectional coupling with *X* as the driver, negative values for *Y* driving *X*, and for symmetric bidirectional coupling.

Since its invention, STE has been used to study interactions in various disciplines ranging from quantum [56], plasma [57] and laser physics [58] via neurology [59], cardiology [60] and anaesthesiology [61–64] to the neurosciences [65]. In addition, extensions have been proposed that allow one to estimate the dominating direction of information flow even from short or transient signals [65]. A recent modification, termed partial symbolic transfer entropy [66,67], extends the STE for multivariate time series and can help distinguish between direct and indirect interactions.

Estimating interaction properties with information-theoretic measures is prone to biases, and statistical errors that depend on the method used and on the characteristics of the data [44,68–71]. For both, the order parameter and the directionality index , asymptotic distributions are unknown, and in the following we apply averaging procedures (both in time and space) to decrease the variance of estimation error and consider both aspects of interaction, strength and direction, to minimize the risk for misinterpretations and to effectively distinguish the various coupling regimes.

Based on previous investigations that employed bivariate time series analysis techniques to characterize *both* strength and direction of interactions in studies on coupled model systems [31,49,50,72,73], we expect the order parameter and the directionality index to depend on the coupling strength *κ* in some driver–responder system as follows (cf. figure 1): increases monotonously with increasing *κ* until for fully synchronized systems. This dependency thus correctly reflects an increasing strength of interaction owing to the increased coupling. On the other hand, correctly detects the direction of interaction for some intermediate values of *κ* only. Moreover, , for uncoupled (*κ*=0 and ), or for fully synchronized systems, (). Distinguishing these coupling regimes thus requires checking the strength of interaction. When investigating interacting systems using empirical data, the coupling strength is usually not known *a priori*, but we expect the aforementioned dependencies to hold (at least approximately) for such data.

## 3. Measuring strength and direction of interactions in the human epileptic brain

We retrospectively analysed strength and direction of interactions between various brain regions by applying the aforementioned methods to multi-channel, multi-day electroencephalographic (EEG) data recorded intracranially in 11 patients (four women, seven men; mean age at onset of epilepsy 10.5 years, range 0–32 years; mean duration of epilepsy 28 years, range 5–53 years). These data had been recorded during the pre-surgical evaluation of intractable focal seizures with different anatomical onset locations. The patients had signed informed consent that their clinical data might be used and published for research purposes. The study protocol had previously been approved by the ethics committee of the University of Bonn. EEG data were recorded from, on average, 48 sites (range 21–73; cf. figure 2), band-pass-filtered between 1 and 45 Hz and sampled at 200 Hz (sampling interval Δ*t*=5 ms) using a 16 bit analogue-to-digital converter. Recordings lasted, on average, 98.6 h (range 25.4–206.2 h) during which seven seizures per patient (range 1–21) were captured. In five patients, seizures originated from the mesial temporal lobe, in two patients from the lateral temporal lobe, and in four patients from the frontal lobe. Post-operatively, all patients were seizure-free [74].

As a compromise between the statistical accuracy for the calculation of interaction indices and approximate stationarity [75], we divided the data into a sequence of non-overlapping segments of 20.48 s duration (corresponding to 4096 data points). This allowed us to calculate and for each combination of pairs of recording sites in a time-resolved manner. If not stated otherwise, all calculations were performed with *m*=5, *l*=3, *N*_{w}=2048 and *N*_{η}=204 [76,77].

Using knowledge concerning location and extent of the SOZ, which is defined by the electrode contacts showing initial seizure activity [10], we assigned all electrode contacts to three location categories:

— focal (): electrode contacts located within the SOZ (on average for all patients 13.0% of all contacts, varying between 5.0% and 20.8%);

— neighbour (): electrode contacts not more than two contacts distant to those from category (on average 9.1%, varying between 2.6% and 14.9%);

— other (): all remaining electrode contacts (on average 77.9%, varying between 66.2% and 87.7%).

This allowed us to concisely analyse interactions between brain regions and to perform inter-subject comparisons, despite the high variability in location and extent of the individual SOZ, which necessitated an electrode implantation tailored to the individual patient.

Because the order parameter is a symmetric and the directionality index an antisymmetric measure, we will present our findings on the strength respective direction of interactions from six combination categories (-, -, -, -, -, -; the direction of interaction is encoded in the sign of ).

In figure 3, we show, as an example, temporal evolutions of the strength () and the direction () of interactions between various brain regions over a time course of 28 h, during which a seizure was captured. The data were calculated from intracranial EEG recordings of the patient whose electrode implantation scheme is shown in figure 2. For only few pairs of recording sites can we observe temporarily increased values of . These increases do not appear to coincide with seizure-related activities, but are nevertheless confined to interactions within the SOZ (-) and to neighbour–neighbour interactions (-). Recognizable indications for directed interactions can be mainly observed for combination categories involving the SOZ, and particularly for interactions between the SOZ and its neighbourhood as well as with other, remote brain regions. These observations indicate the SOZ to be a driving brain region (cf. [30,31,78,79]) with interactions that are not only local but involve large parts of the epileptic network.

From the data shown in figure 3, it can be deduced that, for the most part, indications for directed interactions go along with low-to-intermediate strengths of interactions, whereas high values of the latter are accompanied by indications for symmetric bidirectional couplings. Our findings for all pairs of recording sites (figure 4*a*) and for all combination categories (figure 4*b*) support our expectations mentioned above (figure 1) and underline the importance of investigating both strength and direction of interaction in order to avoid misinterpretations.

We now investigate whether the aforementioned observations of strongest interactions confined to the SOZ and/or its immediate surroundings extend beyond exemplary data. In figure 5, we show and for all categories of combinations of recordings sites estimated from the data from all patients. On this group level, strengths of interactions appear to decrease with an increasing distance from the SOZ and its surroundings, which is in line with a number of previous studies [29,34,80–88]. If we search, however, for the maximum strength of interactions among all combinations of recording sites in each patient separately, we observe this maximum (and the next 10 lower values) for interactions far off the SOZ (category -) in nine of 11 patients. In one patient, categories - and - rank almost equally among the 10 highest values. In another patient, categories - and - clearly rank highest among all combinations of recording sites. These observations on a single patient level are in contrast to many of the previous studies mentioned before, and it remains to be shown whether differences can be related to the much lower number of recording sites taken into account before or to possibly different sensitivities of applied estimators for the strength of interactions [29,41,89].

So far, we considered data only from the seizure-free interval (inter-ictal state) for our analyses. We now proceed by addressing the question whether changes in directed interactions between brain regions constituting the epileptic network can be identified prior to seizures (pre-ictal state). Our previous investigations [90] indicate that directed interactions between the SOZ and homologous contralateral sites decrease (i.e. from a unidirectional driving to a symmetric bidirectional coupling) during the pre-ictal states in 11 of 15 patients with mesial temporal lobe epilepsy. In this study, we also employed symbolic transfer entropy, but restricted our analyses to intrahippocampal recordings. In the following, we will again assume that a pre-ictal phase of 4 h duration exists [21] and compare the distributions of values of from the pre-ictal periods with those from inter-ictal periods (all data that were recorded at least 4 h prior to and 30 min after a seizure) and taking into account all interactions. In figure 6, we show exemplary matrices representing the spatial distribution of averaged strengths and directions of all pair-wise interactions during pre-ictal and inter-ictal periods. For both periods, highest values of can be observed for interactions within the SOZ (category -; right hippocampal formation; sites TR01–TR10) and its surrounding brain areas (category -; sites TBAR1–TBAR3 and TBPR1–TBPR3) but also for interactions in the ipsilateral neocortex, particularly at parieto-occipital sites (category -; TLR06–TLR08 and TLR14–TLR16). Strongest interactions also involve the most anterior contact in the contralateral hippocampal formation (category -; particularly site TL01) and parieto-occipital sites in the contralateral neocortex (category -; sites TLL07–TLL08). In line with Mormann *et al*. [91], we observe—in both hippocampal formations—a gap in the inter-regional strength of interactions namely within the entorhinal cortex and within the hippocampus that indicates the existence of independent rhythms (mostly in the frequency band 0.5–7 Hz) in different subregions of the human medial temporal lobe, produced by autonomous generators. Prominent differences between data from the pre-ictal and the inter-ictal periods are visible for the SOZ, its neighbourhood, homologous sites within the contralateral mesial temporal lobe, and for interactions between inferior and lateral sites from the contralateral temporal cortex. We observe both a pre-ictal increase and decrease of , in line with previous studies that used other, mostly phase-based estimators for the strength of interactions [20,25,92–97].

If we assume that the direction of interactions can best be resolved at intermediate strengths (figure 1), we accordingly observe high absolute values of particularly for interactions involving categories -, - and -. During both the pre-ictal and the inter-ictal periods, the SOZ (right hippocampal formation; sites TR01–TR10) appears to persistently drive all other brain regions, whereas these regions appear to drive structures near the SOZ (sites TBAR1–TBAR3 and TBPR1–TBPR3). Nevertheless, we also observe directed interactions between brain regions that are apparently not associated with the SOZ (category -). As a result, no clear-cut directionality patterns can be observed from the differences between data from the pre-ictal and the inter-ictal periods, although strongest changes in directionality, i.e. more pronounced indications for a unidirectional driving coupling, are confined to brain regions surrounding the SOZ (category -; sites TBAR1–TBAR3 and TBPR1–TBPR3). In figure 7 we show the inter-ictal and pre-ictal frequency distributions of indices and for all categories of combinations of recordings sites. Pronounced relative changes (more than 20%) can be observed only for directed interactions between brain regions surrounding the SOZ and all other sampled regions (category -; inter-ictal: −0.09±0.18; pre-ictal: −0.15±0.16) and for directed interactions between the SOZ and its neighbourhood (category -; inter-ictal: 0.26±0.15; pre-ictal: 0.33±0.13).

In order to investigate whether the aforementioned findings extend beyond exemplary data, we estimated—for each combination category (but excluding those that involve the same location category)—the number of patients for which pre-ictally an increase or a decrease (of at least 5%), or no change of the order parameter and of the directionality index can be observed when compared with the inter-ictal periods (table 1). For about half the cases, the strength of interactions between the SOZ, its surroundings and other brain regions (categories - and -) increased during the pre-ictal periods. The strength of interactions between brain regions not involving the SOZ (category -) remained unchanged pre-ictally in the majority of cases. The direction of interactions increased pre-ictally in about half the cases for interdependencies between the SOZ and its surroundings (category -) as well as between brain regions not involving the SOZ (category -). The direction of interactions decreased pre-ictally in about half the cases for interdependencies between the SOZ and other brain regions (category -).

In table 2, we report the number of patients for which indicates a preferred direction of interactions (preferentially driving versus preferentially responding) during pre-ictal and inter-ictal periods. For this purpose we calculated, for each investigated combination category, the temporal and spatial median of and of (denoted as and , respectively), and regarded the interaction as preferential driving (responding) if was positive (negative), and the accompanying values of were from the interquartile range of the respective distribution of values of , we could neither observe a clear-cut indication for a preferential driving nor for a preferential responding between the investigated brain regions. Pre-ictally, the preferential driving of surrounding brain region by the SOZ (category -; as seen inter-ictally) appears to be lost, however, the driving of surrounding brain regions by other regions (category -) can be observed more often. Although speculative, these observation probably point to an effective *inhibitory surround* [98,99] that prevents the spreading of epileptiform activities during the seizure-free interval.

## 4. Conclusion

Using approaches from symbolic analysis, we investigated—in a time-resolved manner—changes in strength and direction of short- to long-ranged interactions between brain regions constituting the epileptic network. We here employed the order parameter and the directionality index (derived from symbolic transfer entropy), both allowing for a robust and computationally fast quantification of the mentioned interaction properties. We followed previous recommendations [31,50] to investigate both aspects of interactions in order to avoid misinterpretations and to allow an effective differentiation between various coupling regimes.

Analysing multi-day, multi-channel EEG recordings from 11 epilepsy patients, we could confirm some previous findings concerning the strength of interactions that had been derived with frequency-based, phase-based, state-space-based or statistical approaches. Summarizing these observations, the strength of interactions decreases with an increasing distance from the SOZ and its surroundings during the seizure-free (inter-ictal) periods, whereas it exhibits a possibly location-specific increase or decrease during the pre-ictal periods. Our findings obtained from investigating interdependencies between all sampled brain regions, however, indicate that the maximum strength of interactions can be observed for interactions between brain regions far off the SOZ in the majority of cases. It thus remains to be shown whether the dynamics of the SOZ can indeed be characterized by an elevated local synchrony or elevated strength of interactions [80–83,86–88] and whether this property can help to identify the SOZ in the pre-surgical evaluation. It also remains to be shown whether the location-specific changes of the strength of interactions (increase or decrease) during the pre-ictal periods can be regarded as precursors of epileptic seizures. Such an investigation would require using statistical methods for testing the significance of the predictive performance of the applied analysis techniques [25,100–106]. The same holds for the direction of interactions, because the unambiguous inference of the direction of interactions remains an unsolved problem, despite some developments to test the significance of directionality estimates [107–110].

We here considered a spatial average over pre-defined location categories using knowledge concerning location and extent of the SOZ and observed in some patients a preferential driving between some brain regions, whereas in other patients directionality was inverted for the very same brain regions. This ambiguity persisted when comparing data from the inter-ictal and the pre-ictal periods. Owing to these inconsistencies, we cannot yet give a clear-cut answer to the question of who is driving whom in large-scale epileptic brain networks.

There are, by now, only a few studies that investigated directed interaction in the epileptic brain [15,90,111–117], and these studies mostly concentrated on seizures or other epileptiform activities (such as epileptic spikes) and/or were mostly restricted to selected recording sites. A comparison of our findings with those obtained in these studies is thus rather constricted.

The inference of directed interactions from empirical data may be limited by a number of influencing factors. Bivariate analysis, as used here, is often not sufficient to reveal the correct interaction structure, i.e. distinguishing direct and indirect interactions. An extension to the concept of symbolic transfer entropy that aims at accounting for the presence of other observed confounding variables has recently been proposed (so-called partial symbolic transfer entropy [66,67]), but its suitability for characterizing directed interactions in epileptic brain networks remains to be investigated. It is well known that interactions between and within brain regions may be delayed, with time delays reaching up to 200 ms depending on region and functions [118], and not accounting for such delays may lead to spurious indications for directed interactions. An improved characterization of such delayed, directed interactions may be achieved with recently proposed analysis techniques that explicitly incorporate time delays in estimates for information transfer [119–123].

In the future, we expect further improvements in the data-driven characterization of weighted and directed, short- to long-ranged interactions in evolving epileptic brain networks to advance our understanding of the dynamical disease epilepsy, which may guide new developments for individualized diagnosis, treatment and control.

## Funding statement

This work was supported by the Deutsche Forschungsgemeinschaft (grant no. LE 660/5-2).

## Acknowledgements

The authors are grateful to G. Ansmann, S. Porz and C. Geier for critical comments on earlier versions of the manuscript.

## Footnotes

One contribution of 12 to a theme issue ‘Enhancing dynamical signatures of complex systems through symbolic computation’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.