# A Robust Method for Detecting Interdependences: Application to Intracranially Recorded EEG

## 1 Introduction

During the last years the analysis of synchronization phenomena received increasing attention. Such phenomena occur in nearly all sciences, including physics, astrophysics, chemistry, and even economy. Probably the most important applications are in biology and medical sciences. In living systems, synchronization is often essential in normal functioning, while abnormal synchronization can lead to severe disorders. Typical examples are from neurosciences, where synchronization under normal conditions seems to be essential for the binding problem , whereas epilepsies are related to abnormally strong synchronization.

Synchronization can manifest itself in different ways. At one extreme are coupled identical deterministic chaotic systems, which can synchronize perfectly: once the coupling exceeds a critical value, both systems move along identical orbits . If the coupled systems are not identical, in general, they can not move along identical orbits. If they are both chaotic and noise-free, a strict relationship can still exist, provided the coupling is sufficiently strong. Let us denote by $`X=(x_1,\mathrm{},x_N)`$ and $`Y=(y_1,\mathrm{},y_N)`$ two time sequences from which state vectors $`𝐱_n`$ and $`𝐲_n`$ can be reconstructed, e.g., as delay vectors. Let us also assume that one of the systems, say $`X`$, is driving the other. By this we mean that the evolution of $`𝐱_n`$ is autonomous, while $`𝐲_{n+1}`$ is a function of $`𝐲_n,𝐱_n`$, and probably of some external noise . If there is no noise, and if the driving is non-singular, $`𝐲_{n+1}=𝐅(𝐱_n,𝐲_n)`$ with $`det(F_i/x_{nk})0`$, this relationship can always be inverted (at least locally) and can be written as $`𝐱_n=𝐆(𝐲_n,𝐲_{n+1})`$ or, after eventually increasing the embedding dimension of $`Y`$, as $`𝐱_n=\mathrm{\Phi }(𝐲_n)`$ . The opposite relation (probably with some time shift $`k`$)

$$𝐲_n=\mathrm{\Psi }(𝐱_{nk})$$
(1)   
is not guaranteed, although it looks a priori more natural in view of the fact that $`X`$ is assumed to drive $`Y`$. If eq.$`(1)`$ holds for some finite $`k`$, i.e. if the state of the driven system is a unique function of the driver‘s state, this is referred to as ‘generalized synchronization’ . Here two cases have to be distinguished: strong generalized synchronization corresponds to smooth functions $`\mathrm{\Psi }`$, while weak generalized synchronization can lead to functions which may even be nowhere continuous . In the latter case it might be difficult to detect synchronization by observing $`X`$ and $`Y`$, while it is immediately seen when comparing two realizations $`Y^{(a)}`$ and $`Y^{(b)}`$ of the response system: if both are unique functions of the same $`X`$, then obviously $`Y^{(a)}=Y^{(b)}`$, i.e. they synchronize perfectly.

Notice that this notion of ‘generalized synchronization’ is closer to the notion of interdependence, rather than to a mere time shift generating temporal coincidences (this is what the word synchronization actually means). If a softening of the concept of synchronization is accepted in this way, this ‘generalized synchronization’ is clearly not yet the weakest and most general form of synchronization. The weakest form is given just when $`X`$ and $`Y`$, considered as stochastic processes, are not independent. The problem of finding weak effects of synchronization is thus equivalent to find weak interdependences. This is particularly true for a system as complex as e.g., the brain, where the question wether eq.$`(1)`$ holds might be meaningless.

Driver/response asymmetries, as mentioned in the above example, are indeed quite common also in stochastic systems. Distinguishing the driver from the responder is of course one of the central goals, particularly in medicine where it is of utmost importance to detect causal relationships. Unfortunately, no general method exists to detect such relationships unambiguously. Even if $`Y`$ follows the motion of $`X`$ with a time delay as in eq.$`(1)`$, so that $`Y`$ hardly could drive $`X`$, this does not proof that $`X`$ drives $`Y`$. Both systems might be driven by an unobserved third system $`Z`$.

In particular, eq.$`(1)`$ by itself does not imply that $`X`$ drives $`Y`$. This is obvious in cases where $`\mathrm{\Psi }`$ is bijective, i.e. where also $`\mathrm{\Psi }^1()`$ is unique. If $`\mathrm{\Psi }`$ is not bijective (which, as we have seen, actually happens if $`Y`$ drives $`X`$ but fails to synchronize it), then, in general, there are several states of $`X`$ which map onto a single state of $`Y`$. This will typically happen if the state space of $`X`$ is larger than that of $`Y`$. For practical applications where strict equality cannot be observed but only closeness, this means that $`X`$ has a larger attractor dimension (i.e. more effective degrees of freedom) than $`Y`$. But this does not imply any causal relationship.

Typical observables used for detecting interdependences and synchronization are mutual information and cross correlations. Closely related to cross correlations are cross spectra. The main disadvantage of the latter two is that they measure only linear dependences. Causal relationships can (with the above caveats) be tested using time delays, i.e. by comparing $`x_my_n`$ with $`x_ny_m`$. Mutual information is sensitive to all kinds of dependencies (it is zero only if $`X`$ and $`Y`$ are strictly independent), but its estimation imposes quite substantial requirements on the amount and quality of the data. In particular, if the suspected optimal embedding dimension is high, these requirements might be hard to meet. Finally, cross correlations and mutual information are symmetric in $`X`$ and $`Y`$, so that causal relationships can be detected only if they are associated with time delays. A priori, causal relationships might exist without detectable delays and, as we have pointed out, there might exist delays which do not reflect the naively expected causal relationship.

A new class of asymmetric interdependence measures which might overcome some of these limitations has been proposed recently . These authors have assumed that a deterministic relationship as in eq.$`(1)`$ exists, and have therefore not optimized their observables so as to detect reliably weak interdependences in a noisy environment. Moreover, they assumed that eq.$`(1)`$ automatically implies a causal relationship. That this is not unproblematic was discussed above. It is also seen from the fact that the authors of  and  drew exactly the opposite conclusions from mutual predictabilities of $`X`$ and $`Y`$. Equation $`(1)`$ was interpreted in  as indicating that $`Y`$ is the driver and $`X`$ the response, and that $`Y`$ can be better predicted from $`X`$ than vice versa. The opposite interpretation — namely that the response can be better predicted from the driver — was given in . Nevertheless, these observables have been applied successfully to neurophysiological problems .

In the present paper we present another interdependence measure following closely references . But we do not assume eq.$`(1)`$ and we try to make our definition such as to be most robust. Our observable, together with several alternatives, is defined in the next section. Applications to EEG signals recorded from electrodes implanted under the skull of patients suffering from severe epilepsies are presented in Sec.3, while our conclusions are drawn in Sec.4.

## 2 Outline of the Method

Let $`X=(x_1,x_2,\mathrm{},x_N)`$ and $`Y=(y_1,y_2,\mathrm{},y_N)`$ denote two different simultaneously observed time sequences. Typically, they will be measurements of different observables of the same complex system, or measurements taken at different positions of a spatially extended system. The internal dynamics of the system is not known. In particular, it is not known whether the system is deterministic or stochastic, but we are mostly interested in cases where the latter is more likely a priori, or where it is at least unlikely that the attractor dynamics is so low that methods developed specifically for chaotic deterministic systems would be applicable. Physical time is related to the index $`n`$ of $`x_n`$, respectively $`y_n`$ by $`t=t_0+ϵn`$.

Time-delay embedding  in an $`m`$-dimensional phase-space leads to phase-space vectors $`𝐱_n=(x_n,\mathrm{},x_{n(m1)\tau })`$ and $`𝐲_n=(y_n,\mathrm{},y_{n(m1)\tau })`$. The delay $`\tau `$ can be chosen as 1, but for oversampled sequences it might be useful to use some integer $`\tau >1`$. To simplify notation, we assume that also values $`x_{2m},\mathrm{},x_0`$ and $`y_{2m},\mathrm{},y_0`$ are given, so that all delay vectors with index $`1nN`$ can be formed, and the time sequences of delay vectors have $`N`$ elements each. The arrays of all delay vectors will be denoted $`𝐗=(𝐱_1,\mathrm{},𝐱_N)`$ and $`𝐘=(𝐲_1,\mathrm{},𝐲_N)`$.

Let $`r_{n,j}`$ and $`s_{n,j}`$, $`j=1,\mathrm{},k`$ denote the time indices of the $`k`$ nearest neighbours of $`𝐱_n`$ and $`𝐲_n`$, respectively. Thus, the first neighbour distances from $`𝐱_n`$ are $`d(𝐗)_n^{(1)}𝐱_n𝐱_{r_{n,1}}=\mathrm{min}_q𝐱_n𝐱_q`$, $`d(𝐗)_n^{(2)}𝐱_n𝐱_{r_{n,2}}=\mathrm{min}_{qr_{n,1}}𝐱_n𝐱_q`$, etc., where $`𝐱𝐱^{}`$ is the Euclidean distance in delay space, and similar for $`𝐲_n`$. For each $`𝐱_n`$, the squared mean Euclidean distance to its $`k`$ closest neighbours is defined as

$$R_n^{(k)}(𝐗)=\frac{1}{k}\underset{j=1}{\overset{k}{}}\left(𝐱_n𝐱_{r_{n,j}}\right)^2$$
(2)   
while the conditional mean squared Euclidean distance, conditioned on the closest neighbour times in the time series $`𝐘`$, is

$$R_n^{(k)}(𝐗|𝐘)=\frac{1}{k}\underset{j=1}{\overset{k}{}}\left(𝐱_n𝐱_{s_{n,j}}\right)^2.$$
(3)   
Notice that the only difference between these two is that we used the ‘wrong’ time indices for the neighbours in eq.(3). Instead of summing over nearest neighbours, we sum over those points whose equal time partners are nearest neighbours of $`𝐲_n`$. Similarly we define

$$R_n^{(k)}(𝐘)=\frac{1}{k}\underset{j=1}{\overset{k}{}}\left(𝐲_n𝐲_{s_{n,j}}\right)^2$$
(4)   
and

$$R_n^{(k)}(𝐘|𝐗)=\frac{1}{k}\underset{j=1}{\overset{k}{}}\left(𝐲_n𝐲_{r_{n,j}}\right)^2.$$
(5)   
If the point cloud $`\{𝐱_n\}`$ has average squared radius $`R(𝐗)=R^{(N1)}(𝐗)`$ and effective dimension $`D`$ (for a stochastic time series embedded in $`m`$ dimensions, $`D=m`$), then $`R_n^{(k)}(𝐗)/R(𝐗)(k/N)^{2/D}1`$ for $`kN`$. The same is true for $`R_n^{(k)}(𝐗|𝐘)`$ if $`𝐗`$ and $`𝐘`$ are perfectly correlated, i.e. if there is a smooth mapping $`𝐱_n=\mathrm{\Psi }(𝐲_n)`$. On the other hand, if $`𝐗`$ and $`𝐘`$ are completely independent, then $`R_n^{(k)}(𝐗|𝐘)R_n^{(k)}(𝐗)`$. Accordingly, we introduce local and global interdependence measures $`S_n^{(k)}(𝐗|𝐘)`$ and $`S^{(k)}(𝐗|𝐘)`$ as

$$S_n^{(k)}(𝐗|𝐘)\frac{R_n^{(k)}(𝐗)}{R_n^{(k)}(𝐗|𝐘)}$$
(6)   
and

$$S^{(k)}(𝐗|𝐘)\frac{1}{N}\underset{n=1}{\overset{N}{}}S_n^{(k)}(𝐗|𝐘)=\frac{1}{N}\underset{n=1}{\overset{N}{}}\frac{R_n^{(k)}(𝐗)}{R_n^{(k)}(𝐗|𝐘)}.$$
(7)   
Since $`R_n^{(k)}(𝐗|𝐘)R_n^{(k)}(𝐗)`$ by construction, we have

$$0<S^{(k)}(𝐗|𝐘)1.$$
(8)   
If $`S^{(k)}(𝐗|𝐘)(k/N)^{2/D}1`$, then obviously $`𝐗`$ and $`𝐘`$ are independent within the limits of accuracy. If, however, $`S^{(k)}(𝐗|𝐘)(k/N)^{2/D}`$, we say that $`𝐗`$ depends on $`𝐘`$, thereby without implying any causal relationship. This dependence becomes maximal when $`S^{(k)}(𝐗|𝐘)1`$.

The opposite dependences $`S_n^{(k)}(𝐘|𝐗)`$ and $`S^{(k)}(𝐘|𝐗)`$ are defined in complete analogy. They are in general not equal to $`S_n^{(k)}(𝐗|𝐘)`$ and $`S^{(k)}(𝐗|𝐘)`$. Both $`S^{(k)}(𝐗|𝐘)`$ and $`S^{(k)}(𝐘|𝐗)`$ may be of order 1. Therefore $`𝐗`$ can depend on $`𝐘`$, and at the same time can $`𝐘`$ depend on $`𝐗`$. If $`S^{(k)}(𝐗|𝐘)>S^{(k)}(𝐘|𝐗)`$, i.e. if $`𝐗`$ depends more on $`𝐘`$ than vice versa, we say that $`𝐘`$ is more “active” than $`𝐗`$. Again we do not imply this to have any causal meaning, a priori. An important question is whether an active/passive relationship, as defined in this way, has a causal driver/response interpretation in certain circumstances.

In order to understand the origin of active/passive relationships, we consider the simple case where both time sequences are identical, $`X=Y`$, but we use different embedding dimensions $`m_X`$ and $`m_Y`$ in the delay vector construction. More precisely, we take $`m_X<m_Y`$ and $`m_X<m_{\mathrm{opt}}`$, where $`m_{\mathrm{opt}}`$ is an optimal embedding dimension in the sense that for $`m<m_{\mathrm{opt}}`$ the point cloud $`\{𝐱_n\}`$ is not completely unfolded, while it is unfolded for $`mm_{\mathrm{opt}}`$. Thus each $`𝐱_n`$ can be considered as a singular projection of $`𝐲_n`$, $`𝐱_n=\mathrm{\Psi }(𝐲_n)`$ with non-unique inverse $`\mathrm{\Psi }^1`$. Assume now that $`𝐲_s`$ is a close neighbour of $`𝐲_n`$. Then also $`𝐱_s`$ must be a close neighbour of $`𝐱_n`$. But the opposite is not true: Closeness in $`𝐱`$ space does not imply closeness in $`𝐲`$ space. Therefore, conditioning on times $`s`$ where $`𝐲_s`$ are close neighbours of $`𝐲_n`$ has less effect for neighbours of $`𝐱_n`$ than vice versa, and $`S^{(k)}(𝐗|𝐘)>S^{(k)}(𝐘|𝐗)`$. Although this is not a mathematically rigorous argument, it shows clearly that the active/passive relationship, as defined above, mainly reflects the relative number of degrees of freedom and not a driver/response relationship. Systems with many degrees of freedom (high dimensional “attractors”) are more active than those with few.

Notice, however, that $`S^{(k)}`$ is sensitive only to those degrees of freedom which are excited with amplitudes of order $`R^{(k)}`$. The latter depends, among others, on $`k`$ and on $`N`$. The tendency of (weakly) coupled systems to have degrees of freedom which are excited with very small amplitudes is well known . It often leads to wrong estimates of attractor dimensions, and it can make the observable active/passive relationship to depend on parameters such as $`k`$ and $`N`$ . It might be responsible for the contradictory results of .

Before leaving this section, we point out several possible generalizations and alternatives.

(a) Using the same Euclidean distance to define neighbours and in the sums in eqs.(2)-(5) is not necessary. Instead of the geometrical distance, in eqs.(2)-(5) we could have used any other dissimilarity measure between $`𝐱_n`$ resp. $`𝐲_n`$ and the point clouds $`\{𝐱_{r_{n,j}}\}`$ etc.. If we would have used forecasting errors in local forecasts based on these clouds, we would have arrived at interdependence measures very similar to those of . In , also ‘zero time step’ forecasting was studied. This is most closely related to our observables, but it uses only the distance between $`𝐱_n`$ and the center of mass of the point cloud $`\{𝐱_{s_{n,j}},j=1,\mathrm{}k\}`$, while we use all distances $`|𝐱_n𝐱_{s_{n,j}}|`$ individually. It is clear that the latter contains more information, and should therefore be more sensitive.

(b) Instead of using arithmetic averages as in eqs.(2)-(5) and (7), we could have used geometric or harmonic averages. And we could have replaced the average of ratios in eq.(7) by a ratio of (arithmetic, geometric, or harmonic) averages. Again this could severely change sensitivity and robustness. We have not made an exhaustive test of all alternatives, but we checked that the above definitions are more robust than several alternatives. For instance, replacing eq.(7) by

$$S^{(k)}(𝐗|𝐘)^{}\left[\frac{1}{N}\underset{n=1}{\overset{N}{}}\frac{R_n^{(k)}(𝐗|𝐘)}{R_n^{(k)}(𝐗)}\right]^1$$
(9)   
gave much more noisy results in the applications discussed in the next section which were also much harder to interpret physiologically. This is easily understood. In $`S^{}`$, occasional very small values of $`R_n^{(k)}(𝐗)`$ have much more influence than in $`S`$. Such small values are obtained if $`𝐱_n`$ depends abnormally weakly on $`Y`$, which might arise from some perturbation acting at time $`n`$. Thus $`S`$ is more robust against shot noise than $`S^{}`$. We found similar results when using harmonic averages in eqs.(2)-(5). The main difference between the present paper and  is that these authors were interested in the case of noiseless deterministic attractors and strong interdependences where these considerations play no rôle, and they therefore did not try do find the most robust observable. Also, they dicussed only the case $`k=1`$. This gives the strongest signal, but it is also much stronger affected by noise than $`k>1`$. In the following applications we used $`k=10`$ which seemed to give the best signal to noise ratio (see below).

(c) In eq.(6) we essentially compare the $`𝐘`$-conditioned mean squared distances to the mean squared nearest neighbour distances. Instead of this, we could have compared the former to the mean squared distances to random points, $`R_n(𝐗)=(N1)^1_{jn}(𝐱_n𝐱_j)^2`$. Also, let us use the geometrical average in the analogon of eq.(7), and define

$$H^{(k)}(𝐗|𝐘)=\frac{1}{N}\underset{n=1}{\overset{N}{}}\mathrm{log}\frac{R_n(𝐗)}{R_n^{(k)}(𝐗|𝐘)}$$
(10)   
This is zero if $`𝐗`$ and $`𝐘`$ are completely independent, while it is positive if nearness in $`𝐘`$ implies also nearness in $`𝐗`$ for equal time partners. It would be negative if close pairs in $`𝐘`$ correspond mainly to distant pairs in $`𝐗`$. This is very unlikely but not impossible. Therefore, $`H^{(k)}(𝐗|𝐘)=0`$ suggests that $`𝐗`$ and $`𝐘`$ are independent, but does not prove it. This (and the asymmetry under the exchange $`𝐗𝐘`$) is the main difference between $`H^{(k)}(𝐗|𝐘)`$ and mutual information. The latter is strictly positive whenever $`𝐗`$ and $`𝐘`$ are not completely independent. As a consequence, mutual information is quadratic in the correlation $`P(𝐗,𝐘)P(𝐗)P(𝐘)`$ for weak correlations ($`P`$ are here probability distributions), while $`H^{(k)}(𝐗|𝐘)`$ is linear. This might make $`H^{(k)}(𝐗|𝐘)`$ useful in applications.

(d) Instead of eq.(3) we could have defined the time shifted generalization

$$R_n^{(k)}(𝐗|𝐘,l)=\frac{1}{k}\underset{j=1}{\overset{k}{}}\left(𝐱_n𝐱_{s_{n+l,j}}\right)^2,$$
(11)   
with some (positive or negative) integer $`l`$. The idea behind this definition is that it is not clear a priori that $`𝐱_n`$ is most closely related to the simultaneous vector $`𝐲_n`$. Rather, if there are some time delays in generating either $`x_n`$ or $`y_n`$, the ‘natural’ partner of $`𝐱_n`$ might be $`𝐲_{n+l}`$. In this way we can introduce a further element of asymmetry which could give additional hints on causal relationships.

(e) Up to now, we have assumed in general that we use the same embedding for $`𝐗`$ and for $`𝐘`$. This is not necessary, and we could have used a different embedding dimension $`m`$ and a different delay $`\tau `$ for $`𝐘`$. We did not follow this path since $`𝐗`$ and $`𝐘`$ had similar characteristics in examples studied in the next section. But it is worth while to point out that we can use our interdependence measure for pairs of time series with completely different characteristics (amplitudes, spectra, etc.). Dependence does not imply similarity in any sense!

(f) Instead of the Euclidean distance we could have used any other distance in defining neighbourhoods, e.g. the maximum norm.

## 3 Application

### 3.1 Data Acquisition

We analyzed electroencephalographic signals (EEG) that were recorded in patients suffering from pharmacoresistant focal epilepsies. In these patients freedom of seizures can be obtained by resecting the part of the brain responsible for seizure generation. Taking such sort of data is mandatory as part of the presurgical analysis. The sensoring electrodes are left in the brain for typically 2 to 3 weeks. During this time the patients are also watched by video, so that EEG activity can be matched with behavior, and seizures can be identified from either. The analyses reported here were made after surgery had taken place, and after it had become clear from its success whether the localization of the epileptic focus had been correctly predicted.

EEG was recorded from electrodes implanted under the skull, hence close to the epileptic focus and with high signal-to-noise ratio. In particular, we used two types of electrodes: rectangular flexible grids of $`8\times 8`$ contacts placed onto the cortex, and pairs of needle shaped depth electrodes with 10 contacts each, implanted into deeper structures of the brain (see fig. 1).

EEG signals were sampled at 173 Hz using a 12 bit analog-to-digital (A/D) converter and filtered within a frequency band of 0.53 to 40 Hz. The cutoff frequency of the lowpass filter was selected to suppress possible contamination by the power line. For more details on the data and recording techniques, see  and references given therein. The data sets analyzed in this study had a duration of 10 minutes each (cut out from much longer sequences) and were divided into segments of $`T`$ seconds each. Neighbours were searched only within the same segment.

### 3.2 Parameter Selection

As is well known, details of the delay embedding such as choice of embedding dimension $`m`$ and delay $`\tau `$ can be very important. In principle, the theorems of Takens  and Sauer et al.  state that results should not depend on them if data are noiseless and $`N`$ is arbitrarily large, but reality tells different. Many methods have been proposed to find ”optimal” parameter values. However, appropriate choices of $`m`$ and $`\tau `$ strongly depend on specific aspects of the problem at hand (such as noise level, type of noise, intermittency, stationarity, etc.). Thus general recipes which do not take into account these factors can be misleading. This holds true in particular for estimates of $`m`$ based on false nearest neighbours . One of the most popular recipes  for determining the optimal delay $`\tau `$ is based on minimizing the mutual information in a two-dimensional embedding. But in general the same $`\tau `$ does not minimize the mutual information in an embedding $`m`$ $`3`$ dimensions . The same comment applies to estimates of $`\tau `$ from the first zero of the autocorrelation function. Therefore we used none of these a priori estimates of “optimal” embedding parameters in this study. Instead, we approached the problem empirically by calculating $`S^{(k)}(𝐗|𝐘)`$ and $`S^{(k)}(𝐘|𝐗)`$ for different values of $`m`$, $`\tau `$, $`T`$, and $`k`$. In addition, we applied also a Theiler correction  by restricting the nearest neighbour times $`r_{n,j}`$ and $`s_{n,j}`$ to $`|nr_{n,j}|\tau _{\mathrm{Theiler}}`$ and $`|ns_{n,j}|\tau _{\mathrm{Theiler}}`$, and tested several values for $`\tau _{\mathrm{Theiler}}`$. It is of course not feasible to make a systematic search for all possible combinations of these parameters, but we feel sure that our final choices are reasonable and not too far from the optimum. We made these optimizations out of sample, i.e. we used a well understood ‘training’ data set where we could judge the reasonability of our observables by comparing with the medical diagnosis. This training set was not used as test set in any of the subsequent analyses. The “optimal” parameters are $`m=10`$ (embedding dimension), $`\tau =5`$ (delay in units of sampling time), $`k=10`$ (neighborhood size), $`T=10`$ (segment length in seconds), and $`\tau _{\mathrm{Theiler}}=10`$. Indeed, somewhat better results were in some cases obtained with larger $`k`$ (up to $`k=100`$), but we stuck to the above because it was faster without too much loss of significance. The delay $`\tau =5`$ was implemented by simply decimating the time sequences, thereby reducing effectively the sampling rate from 173 Hz to 34.6 Hz. Thus, each segment contained 346 delay vectors.

### 3.3 Data Representation

#### 3.3.1 Depth Electrodes

From the 20 time sequences recorded via the depth electrodes 400 combinations have to be analyzed. Results can be arranged into a $`20\times 20`$ interdependence matrix $`S_{ij}=S^{(k)}(𝐗_i|𝐗_j)`$. We present our results graphically by means of encoding each pixel in a $`20\times 20`$ array using a grey scale. Pixel $`(i,j)`$ is black if $`S_{ij}=1`$ ($`𝐗_i`$ and $`𝐗_j`$ are identical; this happens on the diagonal), while it is white if $`S_{ij}=0`$. The numbering of channels and their arrangement in the matrix are explained in fig. 2.

Quadrants I and IV represent interdependences between signals from the same (left resp. right) hemisphere, while quadrants II and III show interdependences between different hemispheres. More precisely, if a pixel $`(i,j)`$ in quadrant II is darker than its partner $`(j,i)`$ in quadrant III, the region around contact $`i`$ in the right hemisphere is more active than the region around contact $`j`$ in the left hemisphere. Of particular interest are also average values of $`S_{ij}`$, i.e. averaged over a region symmetric under reflection along the diagonal. The average darkness of such a region is a direct measure of its average interdependences with other parts of the brain involved.

A typical example of a grey scale pattern is shown in fig. 3 exhibiting two regions of high interdependence in both the left hemisphere and the right hemisphere. In this case the depth electrodes were not placed in a completely symmetrical fashion. While the electrode in the left hemisphere had 4 contacts in the entorhinal cortex and 6 contacts in the hippocampus, the right electrode had 3 contacts in the entorhinal cortex and 7 in the hippocampus. This difference (confirmed by MRI images) is clearly seen in fig. 3. In addition, there is a stronger interdependence between entorhinal cortex and hippocampus on the left than on the right side, and the left hippocampus can be assumed to be more active than the right one. Interpretations of the latter will be given in sec. 3.4.1.

#### 3.3.2 Grid Electrodes

Since grid electrodes consisted of 64 contacts, it is not very practical to represent the data in the same way as for the depth electrodes. In addition, labeling the contacts by means of a single index will result in a loss of all neighbourhood information, and the patterns would be hard to interpret. A different representation is obtained by displaying each contact as a plaquette of an $`8\times 8`$ matrix, and indicating the activity patterns by arrows connecting these plaquettes . But also such a picture (which is optimal for a small number of electrodes) is too much packed with information for our present applications to be useful.   
We proceeded differently. We first averaged all 60 matrices obtained by cutting the 10 minutes recording into intervals of 10 seconds. The resulting time-averaged interdependences are called $`\overline{S^{(k)}(𝐗_{i_1,i_2}|𝐗_{j_1,j_2})}`$ where $`(i_1,i_2)`$ and $`(j_1,j_2)`$ are the coordinates of the contacts. We next perform a ranking of all entries in the $`64\times 64`$ matrix except the elements on the diagonal. Using the highest one percent of entries after ranking and taking the lower end as a cutoff $`S_c`$, we define for each contact $`(i_1,i_2)`$ an average activity

$$A_{i_1,i_2}=\underset{j_1,j_2}{}\overline{S^{(k)}(𝐗_{j_1,j_2}|𝐗_{i_1,i_2})}\mathrm{\Theta }(\overline{S^{(k)}(𝐗_{j_1,j_2}|𝐗_{i_1,i_2})}S_c)$$
(12)   
and an average passivity

$$P_{i_1,i_2}=\underset{j_1,j_2}{}\overline{S^{(k)}(𝐗_{i_1,i_2}|𝐗_{j_1,j_2})}\mathrm{\Theta }(\overline{S^{(k)}(𝐗_{i_1,i_2}|𝐗_{j_1,j_2})}S_c).$$
(13)   
The cutoff $`S_c`$ is introduced in order to eliminate the effect of contact pairs with very weak interdependence. For these pairs, $`\overline{S^{(k)}(𝐗_i|𝐗_j)}`$ is dominated by noise, and including them would mainly decrease the signal-to-noise ratio.

Using the coordinates $`i_1`$ and $`i_2`$ we can finally represent $`A_{i_1,i_2}`$ and $`P_{i_1,i_2}`$ as $`8\times 8`$ grey scale matrices. Alternatively, we can add them and represent the sum $`A_{i_1,i_2}+P_{i_1,i_2}`$ as a grey scale matrix. An example is given in fig. 4 exhibiting a region with very strong interdependence near the lower right corner. Its interpretation will be given in the next section.

### 3.4 Results

Our results are illustrated by three examples covering lateralization of the focal brain side, precise focus localization in neocortical epilepsies, and changes of interdependences before an impending seizure. These examples are quite typical. A more systematic study involving statistically significant samples is under way and will be presented elsewhere.

#### 3.4.1 First Example

We analyzed 10 minutes of an interictal (seizure-free interval) EEG of a patient suffering from a so called mesial temporal lobe epilepsy. The clinical workup suggested the epileptic focus to be located in the left hemisphere of the brain. We divided the EEG data set into 60 nonoverlapping consecutive 10 seconds segments and calculated a 20 x 20 $`S`$-matrix for each segment as described above. One of these matrices was already shown in fig. 3. This figure is typical for all 60 matrices in showing more interdependences in the left hemisphere than in the right. This concerns both interdependences within the hippocampus, and between hippocampus and adjacent cortex. Indeed, surgery on the left side resulted in complete seizure control of this patient. This suggests that our proposed measure might be able to lateralize the focal side of the brain.

#### 3.4.2 Second Example

We analyzed 10 minutes of interictal EEG data from a patient suffering from a neocortical lesional epilepsy. In this case an $`8\times 8`$ grid electrode was implanted covering the underlying brain lesion. Again the data set was subdivided as in example one.

A typical activity-passivity matrix obtained by means of the procedure described in sec. 3.3.2 is shown in fig. 4. As already pointed out in sec. 3.3.2, we observed highest interdependences in regions near the lower right corner. Indeed, the patient was operated on exactly in this region (which had been identified during presurgical evaluation) and is now free of seizures.

#### 3.4.3 Third example

In contrast to the afore mentioned examples, where we used only EEG recordings from a seizure free interval and averaged the data over time, we now study $`S`$ as a function of time. Our time resolution is again $`T=10`$ sec. Of particular interest are changes of $`S`$ before an impending seizure, as this could finally lead to its prediction .<sup>1</sup><sup>1</sup>1The results of  use a vague definition of the interictal period and might therefore be questionable. But also changes during seizures and during the postictal (after-seizure) period are of interest.

A sequence of interdependence patterns taken before, during and after a seizure is shown in fig. 5. The pattern of interdependences within the right hemisphere remains almost constant, even during the course of the seizure. On the other hand, $`S`$-values of the left hemisphere change dramatically. As confirmed by successful surgery, the left hemisphere was the focal side in this case. During the preictal stage, $`S`$ decreases from a high initial level to almost zero. Notice, that $`S`$ is very low also in quadrant II directly before seizure onset, indicating that the left hemisphere is much less active. In frame #13, shortly before the onset of the seizure, interdependence builds up again on the left side. It reaches its maximum during the seizure and finally declines towards the interictal level.   
This coincides with findings of Lehnertz and Elger  who found reduced complexity before an impending seizure. Notice that “activity” according to our definition essentially depends on the number of excited degrees of freedom, which is exactly what was measured in . The loss of activity before the seizure onset can be interpreted as a more or less hidden pathological synchronization phenomenon. It is assumed that seizure activity will be induced when a ”critical mass” of neurons is progressively involved in closely time-linked high-frequency discharging. This critical mass might be reached if the preceeding level of synchronization decreases, enabling neurons to establish a synchronization which is high enough to finally lead to seizure activity.

At first sight it may therefore seem paradoxical that interdependences decrease before a seizure. But this might indeed be exactly what happens. In a healthy brain a critical mass is never reached because neurons are strongly tied into networks where they communicate with others. A critical stage may be reached when a large population is “idle” and therefore on the one hand uncorrelated with the rest of the network, but on the other hand, easily recruitable for subsequent coherent pathophysiological activity.

## 4 Discussion

We have presented an observable which can detect dependences between simultaneously measured time sequences. It is similar to other synchronization measures proposed recently, but is somewhat simpler and more robust. With the other measures it shares the property of being asymmetric. In principle, it can be assumed that our measure can indicate causal relationships. This might be useful identifying the driver of the two subsystems emitting the sequences. We claim that such information might be obtainable in principle, but the interpretation is subtle and naive arguments can be quite misleading. Nevertheless, this asymmetry is very interesting. It mainly depends on the difference in ‘activity’ which measures the effective number of excited degrees of freedom. This effective number of active degrees of freedom depends on the scales to which the observable is most sensitive. In principle, in an asymmetric driver-response pair the attractor dimension of the response is always at least as high as that of the driver (if both are deterministic), but this might be relevant only at length scales which are too small to be resolved practically.

Our measure could also be used to detect generalized synchronization, but we do not assume in our applications that the signals are chaotic with low dimensions. In contrast to recent attempts to detect phase synchronization in brain signals, our measure does not treat phase information different from amplitude information, and thus we cannot discuss phase or frequency locking.

We applied our measure to intracranial multichannel EEG recordings taken from patients suffering from severe epilepsies. We found significant dependences between different recording sites, and these dependences were in general not symmetric. Due to the careful pre-operational screening of these patients and their observation after being operated, we could compare our results in detail with other neurophysiological findings. The most interesting preliminary results are the following:

1) During seizure-free intervals, the seizure generating area of the brain exhibited higher interdependences than other brain areas.

2) Some seizures analyzed here were preceeded by short periods (30 s to several minutes) during which extremely low dependences were confined to the seizure generating area.

Although these results are very encouraging, a more systematic study is needed and is under way. In addition, a host of further investigations is imaginable. Obvious candidates are the influences of drugs, the effect of mental activity (epilepsy patients behave normal even with implanted electrodes), or of various stimuli. Another important problem is the determination of the ‘critical mass of neurons needed to trigger a seizure. Moreover, a more systematic comparison with other diagnostical tools is necessary beforehand. Finally, the present findings already suggest a number of physiological results whose interpretation demands a thorough theoretical study. For instance, it is a priori not clear whether a seizure is primarily triggered by a change of activity in the seizure generating area, or a change of susceptibility of the surrounding regions. We hope that the near future will show progress along these lines.

Acknowledgements

We thank J. Müller-Gerking, R. Quian Quiroga, T. Schreiber, and W. Burr for the inspiring discussions and helpful comments during the study.

Figure captions:

Fig. 1: Schematic view of the two types of intracranial electrodes used in this paper. Grids were placed onto the cortex and have either $`8\times 8`$ electrodes. Needle shaped depth electrodes have ten contacts each and were always used pairwise in a left-right symmetrical fashion. In some cases, depth electrodes and grids were used together.

Fig. 2: Scheme of subdivision of the 20x20 matrix $`S_{ij}`$. The indices $`L_1`$ to $`L_{10}`$ denote the contacts on the left depth electrode, from innermost ($`L_1`$) to outermost ($`L_{10}`$). Similarly, $`R_1`$ to $`R_{10}`$ correspond to the right depth electrode. The index $`i`$ runs horizontally, while $`j`$ runs vertically. E.g., quadrant II shows the effect of conditioning right hemispheric channels on the channels from left hemisphere.

Fig. 3: Example for a 20x20 $`S`$-matrix of a 10 sec. segment recorded during the seizure free interval using 10 depth electrodes on each side of the brain.

Fig. 4: (A) Average activity pattern in an $`8\times 8`$ grid electrode; (B) average passivity and (C) normalized sum of both.

Fig. 5: Sequence of interdependence patterns $`S_{ij}`$ including preictal (1-14), ictal (15-16) and postictal (17-20) brain electrical activity.

Figure 1   
Arnhold et al.

Figure 2   
Arnhold et al.

Figure 3   
Arnhold et al.

Figure 4   
Arnhold et al.

Figure 5   
Arnhold et al.