## Abstract

For the past decade, the detection and quantification of interactions within and between physiological networks has become a priority-in-common between the fields of biomedicine and computer science. Prominent examples are the interaction analysis of brain networks and of the cardiovascular–respiratory system. The aim of the study is to show how and to what extent results from time-variant partial directed coherence analysis are influenced by some basic estimator and data parameters. The impacts of the Kalman filter settings, the order of the autoregressive (AR) model, signal-to-noise ratios, filter procedures and volume conduction were investigated. These systematic investigations are based on data derived from simulated connectivity networks and were performed using a Kalman filter approach for the estimation of the time-variant multivariate AR model. Additionally, the influence of electrooculogram artefact rejection on the significance and dynamics of interactions in 29 channel electroencephalography recordings, derived from a photic driving experiment, is demonstrated. For artefact rejection, independent component analysis was used. The study provides rules to correctly apply particular methods that will aid users to achieve more reliable interpretations of the results.

## 1. Introduction

For the past decade, the detection and quantification of interactions within and between physiological networks has become a priority-in-common between the fields of biomedicine and computer science. Prominent examples are the interaction analysis of brain networks [1] and of the cardiovascular–respiratory system [2–5]. In the brain, local neuronal structures fulfil functions as a single system (functional segregation) but are just as important as part of an interconnected, distributed, hierarchical set of network clusters. Interactions in and between such distributed network clusters enable the functional integration and provide, together with the anatomical connectivity, the phenomenal functionality and efficiency of the brain. Current approaches to interaction analysis are linked to the concepts of functional and effective connectivity.

Functional connectivity describes the temporal correlations between neurophysiological events of spatially neighbouring or remote neuronal structures. Effective connectivity is defined as the directed influence of one neuronal structure on another, mediated directly or indirectly. The more prominent concepts of effective connectivity analysis include: dynamic causal modelling [6], transfer entropy [7], Granger causality (GC) [8], directed transfer function (DTF) [9] and partial directed coherence (PDC) [10] which all use different types of functional measurements as input, in particular, electroencephalography (EEG), magnetoencephalography (MEG) and functional magnetic resonance imaging (fMRI). Brain interactions are commonly transient, and several groups of neurons interact with each other either sequentially and/or simultaneously. Therefore, *time-variant multivariate methods* are required for the identification, quantification and modelling of the spatio-temporal dynamics of such complex interaction networks. The time-variant Granger causality index (tvGCI) and PDC approaches on the basis of a time-variant, multivariate autoregressive (AR) model have been introduced, evaluated and applied by authors working together from Jena and Rome. The applications are related to sensors (electrodes) [11] and source space signals by using EEG [12] and BOLD (fMRI) [13] as well as signals derived from the cardiovascular–cardiorespiratory system [2–5,14]. In most EEG/MEG applications, the frequency information is indispensable, i.e. it is of particular interest to investigate superimposed oscillatory brain activities [15,16]. Therefore, time-variant PDC (tvPDC) and PDC-based measures are preferable because of their frequency resolution. Thus, results of connectivity analyses encompass spatial information (multivariate methods) and information with regard to the direction of the interaction in the time as well as frequency domains. Both PDC and the DTF are insensitive to zero-phase delay between two signals [17] which occurs as one effect of volume conduction when connectivity analysis is applied to EEG signals. However, volume conduction results also in a weighted mixing of different source signals, which, in turn, directly affects the estimation of the AR parameters and with it the transfer function which is used for PDC and DTF approaches. This was demonstrated by Gomez-Herrero *et al.* [18], who modified the analysis procedure by a de-mixing of the EEG activity.

We developed two algorithms to estimate the coefficients of time-variant multivariate AR models for GC and PDC analyses, a recursive least-squares (RLS) algorithm [19] and a Kalman filter approach [20]. Both are linear AR model approaches (a review is given in [21]). With the introduction of the new methods, the testing for reliability and a comparison of performance with state-of-the-art methods must be carried out. Thus, specific factors that influence connectivity results were investigated: (i) the estimation method for the time-variant AR coefficients (RLS versus Kalman filter) [20], (ii) the intrinsic and determined model-related parameters (dimension, model order, adaptation parameters for the RLS and the Kalman filter approach [20]), and (iii) the measurement conditions as well as signal characteristics and forgetting factor for the RLS method [12].

Comparable investigations were carried out on time-invariant approaches [22] where the length of the analysis interval must be additionally taken into consideration as an influencing parameter. Time-invariant versions of different connectivity measures have demonstrated reliability in detecting neural connectivity [23,24]. Furthermore, the impact of large differences between the variances of the involved signals was demonstrated by Baccala & Sameshima [10], and generalized PDC was introduced to minimize this effect [25]. In addition, the influence of signal filtering (broadband and notch filtering for artefact rejection) on GC has been thoroughly investigated by Florin *et al.* [26] as well as by Barnett & Seth [18]. Filtering is a frequently used pre-processing step in the analysis of directed interactions. De-noising is another pre-processing procedure; its influence on neural signals was investigated by Nalatore *et al.* [27,28]. The influence of volume conduction on the reliability of the transfer entropy-based interaction analysis was studied by Lindner *et al.* [29].

Independent component analysis (ICA) is frequently applied to reject artefacts from EEG/MEG which originate from non-brain sources (interferences from eye movements, the cardiac field, myoelectric activity, magnetic resonance (MR) gradient switching) [30].

All these factors influence to a greater or lesser extent the results of connectivity analyses when using time-invariant as well as time-variant approaches. Therefore, it is important to define the general requirements for their reliable application and to provide information on limitations and possible pitfalls to the method's user to avoid wrong or misleading interpretations of the results.

This study is a systematic and comprehensive investigation of the impact of various influencing factors on tvPDC analyses by using the multi-trial Kalman filter approach (general linear Kalman filter algorithm). For this, signals were used which encompass an abrupt switchover from one stationary connectivity structure to another. By means of the switchover, the dynamics of the time-variant estimation approach can be analysed. This simulation strategy was adapted to EEG data which are derived from our previous studies investigating the methods [31–33] and used to study entrainment effects on brain activity [34]. The EEG data are from an experiment in which repetitive flicker stimulation was applied (photic driving experiment [34]). In this experiment, periods of time with or without a flicker stimulus alternated (every 4 s), and thus it could be expected that after the onset or end of the stimulation reorganizations of connectivity networks would occur to a certain extent.

The study is organized as follows: in §2, the mathematical fundamentals of the time-variant multivariate AR model are given, and the tvPDC computation is introduced. Additionally, the simulated connectivity networks and EEG data are described. In the following section, the influences of basic estimator and data parameters are demonstrated. Furthermore, the effects of band pass filtering and volume conduction are shown. These investigations were all performed on simulated signals. The results of PDC analysis of the EEG are described in §3*b*. Thereafter, analysis both with and without electrooculogram (EOG) artefact rejection changes was carried out. In this way, the influence of an ICA-based artefact rejection can be exemplarily illustrated. The application conditions are summarized, and the meaning of our results for other methodological approaches and signals is discussed in §4.

## 2. Methods and material

### (a) Methods

#### (i) Time-variant partial directed coherence

For the computation of the tvPDC, a time-variant multivariate AR (tvMVAR) model was used. For *n*=1,…,*N*, let
2.1be a *D*-dimensional tvMVAR process of order *p* with **Y**(*n*)∈R^{D}, **A**_{r}(*n*)∈R^{D×D} and **E** an uncorrelated *D*-dimensional Gaussian process with zero mean. The notion of PDC is usually defined by means of the inverse of the transfer function of the process **Y**. Let
2.2be the Fourier transform of the time-variant AR parameters. Assuming a covariance matrix equal to a multiple of the unit matrix, the tvPDC is defined by
2.3with the complex conjugate .

#### (ii) Filtering

Filtering operations are commonly used as pre-processing steps. To investigate filter effects, we applied different filter types with different band widths to the simulated data. We used acausal FFT-, FIR- and IIR-filters with band limits [0.10 0.15], [0.05 0.20] and [0 0.25] each (normalized frequencies, where *f*=0.5 corresponds to the Nyquist frequency). The middle frequency of all filters equals *f*=0.125, which represents the data peak frequency (figure 1) and the relevant frequency of the PDC analysis. The FIR-filter order was 50 for both band pass filters, and 30 for the low pass filter. All IIR-filters were Butterworth filters of order 3.

#### (iii) Rejection of eye movement artefacts from the electroencephalogram

The elimination of artefactual interferences caused by electrical brain-independent sources from brain activity (e.g. blinking, eye movement, heart beat and muscle activity) is another important pre-processing step in EEG analysis. ICA is a blind source separation method which has become the preferred tool for artefact elimination. One reason is that the superposition of electrical signals via volume conduction fits the linearity assumption of ICA with regard to the mixing process very well. The measured data were decomposed into statistically independent components (ICs), where it is assumed that some of the components only contain non-brain activity. These components have to be identified and removed before a complete projection of the remaining components (EEG activity) is carried out which then yields an artefact-free EEG.

In this study, the algorithm of choice was fastICA, because it is generally agreed that fixed-point algorithms, i.e. fastICA and its derivatives, deliver a good performance on average [35,36].

Various authors favour including the EOG as a measured signal for an improved identification of the brain-independent portion of the EEG. We considered both alternatives to investigate the effects of an EOG inclusion, namely the EEG data are

— fed into an artefact rejection processing based on an ICA (ICA

^{−}) and— fed together with the EOG into an artefact rejection processing based on an ICA (ICA

^{+}).

Additionally, we considered an analogous processing schema without ICA (ICA^{Ø}; figure 2). ICA artefact rejection was applied to the entire EEG recording (20 stimulus trains and the corresponding resting intervals, see §2*b*(ii)) of a defined stimulation frequency. The selection of the IC that most likely corresponds to a non-brain origin may be performed either manually or automatically. We used the latter approach based on maximizing the absolute correlation with the EOG.

#### (iv) Statistical testing

Identification of significant directed interactions was achieved using a bootstrap approach following the concept proposed by Sato *et al.* [37]. Sato *et al*. introduced a non-causality property from component **Y**_{j} to component **Y**_{i} by setting all AR parameters *A*_{r,ij}(*n*)=0. Based on estimated and modified AR parameters as well as resampled residuals, a multitude of PDCs may be realized under the null hypothesis H_{0} of non-causality. Finally, the distribution of *Π*_{i→j}(*f*, *n*) under H_{0} may be estimated based on the bootstrap sample . In the time-variant case, the residual resampling has to be modified, because the distributions of **E**(*n*_{1}) and **E**(*n*_{2}) might be different for *n*_{1}≠*n*_{2}. Thus, a residual mixing with respect to the time *n* is no longer allowed. However, in the presence of multiple trials **y**_{k}, 1≤*k*≤*K*, of a process **Y**, the stationary approach [37] may be generalized to the time-variant situation. In fact, multiple trials provide *K* estimated residual vectors **e**_{k}(*n*), where the corresponding random vectors **E**_{1}(*n*), … ,**E**_{K}(*n*) are independent and identically distributed for each *n*. Consequently, the residual resampling may be carried out with respect to the trial index *k* rather than with respect to the time index *n*.

The agreement of two networks was quantified by Cohen's kappa *κ* [38]. Here, each network is considered as a directed graph and described by a *D*×*D* matrix containing zeros (no directed interaction/edge) and ones (directed interaction/edge). The matrix diagonal is completely equal to zero, because self-loops are not considered. Consequently, the agreement of *D*⋅(*D*−1) binary variables is quantified by *κ*.

A common significance level of *α*=0.05 was used for all statistical investigations. Alpha was adjusted for multiple hypotheses testing by means of the Bonferroni correction.

### (b) Materials

#### (i) Simulated dataset

To demonstrate the effect of different factors on a tvPDC connectivity analysis, we used a consistent simulated dataset. For it, we temporally concatenated two five-dimensional network structures introduced in Baccala & Sameshima [10]. Following our EEG data, the signals were simulated with 20 trials (*K*=20) and 800 samples (*N*=800), where the first 400 samples correspond to the AR process described in fig. 2 of Baccala & Sameshima [10], and the last 400 samples were simulated according to the AR parameter defined in fig. 4 of Baccala & Sameshima [10]. The switch from one network structure to the other and the corresponding amplitude spectra are shown in figure 1. All analyses based on the simulated datasets were performed using 50 realizations (with 20 trials, five channels and 800 samples each).

#### (ii) Electroencephalography datasets

In order to investigate effects of ICA pre-processing, we used a real EEG dataset derived from a photic driving experiment. It is described in detail by Schwab *et al.* [34]. For this photic driving experiment, 10 healthy volunteers (aged 22–40 years, five male, five female) were exposed, eyes closed, to stimulation by an intermittent flickering light while EEG (30 EEG channels, enhanced 10–20 system with a 10–10 system over the occipital region (figure 3*d*), plus one channel EOG, plus one channel EMG (Compumedics Neuroscan, El Paso, TX)) and MEG (31 channels, Philips, Hamburg, Germany) were recorded simultaneously. Data were sampled at 1000 Hz, and hardware filtered between 0.1 and 300 Hz. The data were downsampled to 200 Hz and re-referenced to the common average reference for connectivity analysis as recommended for steady-state visual evoked potentials by Srinivasan *et al.* [39].

The stimuli were delivered via optical fibres from two light-emitting diodes, 9 cm away from the closed eyes of the subjects. Resting periods (30–60 s) were recorded between the stimuli. Each stimulation frequency was presented in a sequence of 20 trains (*K*=20 trials). A single train contained 40 flashes and was separated again by a resting period (see [34] for further details).

In this study, we present the results for one volunteer (no. 6) [34] during repetitive flicker stimulation (individual alpha rhythm frequency: 10.8 Hz). A 11.34 Hz flicker frequency was used. The signal associated with the electrode T8 had to be excluded because it was corrupted by excessive noise. Thus, 29 EEG channels and one single EOG channel were finally processed.

For tvPDC connectivity analyses, we used data segments of 4 s duration: each segment was triggered to the flicker onset, and covers 2 s pre- as well as 2 s post-flicker onset. As we were primarily interested in alpha-band frequencies, we evaluated data only in the 8–13 Hz range. For a consolidated analysis, we aggregated PDC values of the alpha-band by taking the median *π*_{i→j}(*n*)=median_{8≤f≤13} {*π*_{i→j}(*f*,*n*)}.

In order to assess significant interactions, we used the bootstrap method described in §2*a*(iv) with 1000 bootstrap replications.

## 3. Results

### (a) Influence of Kalman filter settings, autoregressive model order, signal-to-noise ratio, filtering procedures and volume conduction—simulated dataset

#### (i) Kalman filter settings

The Kalman filter [20] involves two control parameters. One parameter (*c*_{1}) controls the adaptation of the covariance matrix, where low values imply a smoother temporal course and a slowed reaction to rapid changes of the data's covariance structure. The second parameter (*c*_{2}) defines the step-width of the random walk that is used to update AR parameters. Analogous to *c*_{1}, lower values imply a smoother temporal course of AR parameters at the expense of an increased adaptation time to rapid changes of process characteristics. The parameters should be properly chosen to achieve a compromise between fast adaptation and smoothness of the estimated AR parameters. Figure 4*a* shows the mean-squared model residuals (MSMRs) for different settings of *c*_{1} and *c*_{2}. Obviously, MSMRs (as a function of *c*_{1} and *c*_{2}) develop a flat minimum in the present example. Here, it is remarkable that the structure (in particular, the location of the flat minimum) of MSMRs was very similar to the mean-squared deviation of the estimated and true AR parameters, which actually represents the intrinsic effect measure. However, this measure is not available for real-world data, and thus MSMRs have to serve as an appropriate substitute. Figure 4*b* shows the mean-squared deviation of the estimated and true AR parameters for different settings of *c*_{1} and *c*_{2}. The MSMRs and mean-squared AR parameter deviations do not provide identical but contiguous minima. Thus, owing to the flat minimum, MSMRs serve as a suitable criterion to find adequate values for *c*_{1} and *c*_{2}.

All subsequent analyses of the simulated datasets were performed with *c*_{1}=0.01 and *c*_{2}=0.02.

#### (ii) Autoregressive model order

In the stationary case, the influence of the AR order is well known [38]. In summary, a model order which is too low must be avoided as the AR process properties cannot be correctly considered. A model order which is too high has a negative influence on the estimation quality, and common significances are lost because of too many AR parameters that have to be estimated. Between both extremes, there is usually a broader range, where different model orders result in almost identically identified networks (e.g. as shown in Babiloni *et al*. [40]). We could substantiate these findings in the time-variant case (figure 5*a*).

#### (iii) Signal-to-noise ratio

As shown in Astolfi *et al.* [12] and Hemmelmann *et al*. [13], connectivity analyses based on tvMVAR models are very sensitive with respect to additive noise. To investigate the noise influence on our dataset, we examined different signal-to-noise ratios (SNRs) defined by
3.1where *σ*^{2} denotes the variance of the added zero mean Gaussian noise. In our dataset, the numerator of equation (3.1) was equal to 4.89 and represents the average power with respect to all channels *y*_{1}, … ,*y*_{5}. The SNR represents a global SNR, where the channel conditional SNRs vary owing to different channel powers (see also figure 1). As shown in figure 5*c*, even a slight noise level reduces the network detection accuracy.

#### (iv) Filter pre-processing

In combination with GC measures, filtering operations can significantly influence results of connectivity analyses [26]. Similar effects can be expected for PDC-based analyses, because filters may substantially affect estimated AR parameters. In our dataset, all filters result in a significant impairment of network identification, where the deterioration correlates with the band width and may be observed for all filter types used. As a consequence, it is advisable to avoid narrow-band filtering procedures before estimating the AR parameters.

#### (v) Aspects related to volume conduction

A weighted spreading of one source to two electrodes (sensors S1 and S2) *s*_{1}(*n*)=(1−*ε*)⋅*y*(*n*) and *s*_{2}(*n*)=*ε*⋅*y*(*n*) (0<*ε*<1) basically violates the presupposition of uncorrelated residuals. Nevertheless, the Kalman filter algorithm was able to fit a time-variant AR model, regardless of which *ε* has been chosen. Furthermore, PDC values were estimated correctly, i.e. no coupling between *s*_{1} and *s*_{2} was indicated.

To investigate the relationship between linear mixing and interaction detection, we considered a subnet of our artificial network (figure 1) consisting only of node 1 and node 5. It was used as a model for source activity. A weighted mixing of two coupled sources to two sensors S1 and S2 (figure 6) *s*_{1}(*n*)=(1−*ε*)⋅*y*_{1}(*n*)+*ε*⋅*y*_{5}(*n*) and *s*_{2}(*n*)=(1−*ε*)⋅*y*_{5}(*n*)+*ε*⋅*y*_{1}(*n*) (0<*ε*<0.5) results in two correlated components. When increasing *ε*, the number of false-positive detections increased for both couplings 1→5 (figure 6*a*,*b*) and 5→1 (figure 6*c*). On the other hand, the detection of the existing interaction 5→1 from sample 401 to sample 800 was not negatively influenced by the mixing procedure (figure 6*d*). It was shown that the detection quality can be improved by a prepended principal component analysis (PCA) applied to sensor signals [18,41]. Here, we used PCA with complete variance explanation. The connectivity analysis was carried out in component space, where components and sensor data were associated by means of a maximum correlation between them. After PCA, the number of false-positive detections associated with the non-existent coupling 1→5 is significantly reduced (figure 6*e*,*f*), whereas detection of the existing interaction 5→1 (samples 401–800) remains unchanged (figure 6*h*). By contrast, the false-positive rates with respect to 5→1 (samples 1–400) are impaired even for smaller *ε* (figure 6*g*). In summary, in our example, volume conduction results in a significant increase in false-positive detections for mixed sources. False-negative detections play a secondary role here, which is most likely due to a large data amount regarding time-series length and number of trials.

The proportion of false positives depends mainly on the entire network connectivity structure (including coupling strengths) and the mixing parameter *ε*. An *ε*>0 causes correlated residual noise processes, which violate a fundamental assumption for (stationary) AR processes. Thus, it is not surprising that an AR-based connectivity analysis fails when *ε* becomes too big. The transition from sufficient performance to full failure occurs abruptly, where the change depends on various factors such as the intrinsic coupling situation, the estimator and available samples. The implications for practice are that it is necessary to analyse correlations between single sensor channels or to apply tests for instantaneous mixing before a (AR-based) connectivity analysis is initiated. If needed, alternative approaches should be used [29].

### (b) The influence of eye movement rejection on real electroencephalography data connectivity analysis

The influence of EOG artefacts strongly depends on the particular recording situation and conditions, e.g. analysis of ongoing EEG in healthy subjects or patients, event-related or -evoked potentials with different stimulation paradigms. Therefore, a variety of detection and rejection approaches exists. Our aim is to sensitize the user to the fact that problems in the interpretation of the results may occur when artefact rejection algorithms have been applied. Therefore, the ICA-based rejection of the EOG artefacts should exemplarily demonstrate that a careful handling of such pre-processing procedures in the framework of connectivity analysis is obligatory.

We defined the IC with the maximal absolute sample correlation value as an EOG artefact and used it for correction further on. The absolute correlation value |*r*_{i,EOG}| and the identified component *s*_{i} itself depend on whether the EOG was used for ICA or not (|*r*_{i,EOG}|=0.815 for ICA^{−}, |*r*_{i,EOG}|=0.931 for ICA^{+}).

Following §3*a*(i), the Kalman filter parameters were set to *c*_{1}=0.005 and *c*_{2}=0.01, regardless of which EOG pre-processing variant (ICA^{ø}, ICA^{+}, ICA^{−}) was selected. The tvMVAR order *p* was chosen according to the Akaike information criterion (AIC), and it was tuned to achieve a coincidence between the estimated parameter spectrum and the Fourier spectrum of the signal. The tuning procedure considers that the AIC often results in an overestimated model order. Finally, we used a common model order of *p*=16.

The comparison of ICA^{ø}, ICA^{−} and ICA^{+} shows that the time courses of corresponding PDC values are in most cases very similar. This applies, in particular, in the presence of very high PDC values (figure 3*b*). Nevertheless, there were also a few electrode combinations where we found considerable differences in the PDC values (figure 3*a*,*c*). They are mainly located in the frontopolar and left parietal/occipital regions (grey fields in figure 3*d*).

For the sake of clarity, in addition to the full time resolution of *π*_{i→j}(*n*), we defined averaged PDC values within eight time intervals *I*_{1}, … ,*I*_{8} of 500 ms each. In the bottom line of the processing scheme in figure 2, the adjacency matrices of significant PDC values in time interval *I*_{8} are presented. It becomes apparent that the rough network topology is quite similar for all ICA^{ø}/ICA^{−} and ICA^{+} approaches. A topological view of the significant PDC values is given in figure 3*d*. It shows the electrode scheme with significant PDC values, which are displayed by arrows between the respective channels. This is the scheme for the version without ICA (ICA^{ø}), which mostly agrees with the other two conditions (ICA^{−} and ICA^{+}).

The comparisons ICA^{ø} versus ICA^{−}, ICA^{ø} versus ICA^{+} and ICA^{−} versus ICA^{+} by Cohen's kappa demonstrate that accordance between the different results only slightly varies with time and that the best agreement exists between ICA^{ø} and ICA^{−}, where the largest deviations appear between ICA^{−} versus ICA^{+}.

Temporal trends arising from the photic driving were minimal, even though there were a few PDC time courses which appear to fit to the experimental design (i.e. increase or decrease along with flicker onset). This held true for all three EOG rejection approaches. An example for such seemingly flicker-related PDC courses from ICA^{ø} is shown in figure 7*a*. The first three decline during photic driving, whereas the other three rise. Observing the topographical representation of the depicted channel combinations (figure 7*b*), it is notable that all six originate from the parietal/occipital regions, which makes sense as these regions are known to be involved in visual processing. However, such interpretations remain speculations, because pure PDC values are not meaningful; valid statements require networks based on statistical testing.

## 4. Discussion

A thorough evaluation of a new method by a developer is only one basic requirement in obtaining reliable results on application of the particular method to real-world data. It is also important to outline the requirements for the appropriate application of the method to real-world data. Our study describes limitations in applicability of time-variant connectivity analyses and provides recommendations for their application.

The proper parameter settings (*c*_{1}, *c*_{2}) of the Kalman filter yield a suitable compromise between fast adaptation and smoothness of estimated AR parameters. The MSMR for different settings of *c*_{1} and *c*_{2} provides an appropriate substitute for the lacking intrinsic effect measure.

The well-known rules for the model order selection of stationary AR processes proved to be transferable to the time-variant case: a suitable choice has to balance between appropriate reproduction of data and the estimation quality. A high number of trials is important to improve the estimation accuracy and to increase the adaptation speed (figure 5*b*) [12]. It should be noted that our approaches do not use averaged EEG/ERP/EP data to improve the SNR, because the multi-trial Kalman filter approach considers each realization of the underlying stochastic process separately without prior averaging.

The SNR is a crucial factor which may strongly influence the reliability of the approach. In our study, SNRs less than 10 dB yield an inferior agreement between identified and ground truth networks. The SNR of the data used should be estimated before application. A de-noising of the data improves the reliability of the method as shown by previous studies [27,28].

In the case of tvGCI, it may be desirable to apply a filtering to investigate connectivity in specific frequency bands. However, our results show that the application of filters may significantly affect connectivity patterns and should thus be avoided as pre-processing steps for GC analysis [18,26]. However, such band filtering is not required for PDC calculation, because it is defined in the frequency domain. A data-driven signal decomposition (e.g. empirical mode decomposition) with a subsequent PDC estimation was successfully applied for revealing connectivity in different frequency bands [42]. In comparison with filtering strategies, data-driven signal decompositions might be more efficient pre-processing approaches. From this point of view, PDC is preferable to GCI, if one is particularly interested in the frequency properties of connectivity patterns.

Volume conduction induces correlated sensor activities at neighbouring sensors by superposition of underlying brain source activities. In order to reduce these effects, Marzetti *et al.* [43] used source PCA as a pre-processing step to separate uncorrelated sources. Subsequently, pairs of correlated sources were unmixed by minimum overlap component analysis [44]. Gomez-Herrero *et al.* [18] used PCA and source modelling in parallel as pre-processing steps.

Analyses at sensor level suffer from a limited capability to interpret the identified connectivity patterns. In particular, this is due to volume conduction effects which may mask or pretend interaction between neuronal sources at sensor level. Additionally, identified interactions between sensor signals cannot be explicitly mapped to interactions between cortical sources. On the other hand, the application of connectivity measures to source activities is based on the adequateness of specific modelling assumptions (source and head model). Ultimately, analyses in the source space seem to be the most promising approaches because of the aforementioned limitations of sensor-based analyses.

An ICA-based EOG artefact removal may change the number of significant connections in comparison with connectivity networks computed without EOG rejection. The restriction of the effects to a few electrodes is most probably attributed to the frequency range used in our study which is 8–13 Hz, whereas ocular artefacts show a peak at 2 Hz [45]. Additionally, the gamma frequency range can be contaminated by artefacts which originate miniature saccades [46]. Such induced high-frequency artefacts can be identified by tvPDC analysis between EOG and EEG/MEG, as previously demonstrated in Wacker *et al*. [31], and rejected by an ICA approach [47]. Artefacts incorporating effects of moving sources as, for example, heart-related artefacts cannot be fully tackled by ICA. Thus, pure ICA-based approaches for heart beat artefact removal might miss the moving artefact part, which may result in spurious interactions at sensor level owing to varying arrival times of heart-beat-related signal components at different sensors. Here, ICA in combination with a subsequent template-based artefact removal might be an appropriate approach.

A linear approach cannot identify nonlinear interactions. However, linear AR models can be extended by nonlinear terms. A generalized GCI based on multivariate self-exciting threshold AR (SETAR) models was introduced in Leistritz *et al*. [48]. The SETAR approach is sensitive to a wide range of nonlinear signal properties. Non-AR-based approaches were developed to detect nonlinear interactions [49]. Transfer entropy [2,7,50] represents a kind of nonlinear GC, measuring the deviation from an assumed independency of two systems. Equivalence between GC and transfer entropy for Gaussian variables was recently demonstrated [51]. A multivariate version of transfer entropy was introduced [2,50].

All the influencing factors discussed above also affect the estimation quality of nonlinear non-AR approaches. Yet, many of the nonlinear methods are restricted to bivariate processes and require longer analysis intervals than AR-based methods [52]. Consequently, it is important to know the advantages and disadvantages of the methods in relationship to the properties of the signals which are used to analyse connectivity.

## Acknowledgements

The study was supported by DFG grant no. Wi 1166/9-2.

## Footnotes

One contribution of 13 to a Theme Issue ‘Assessing causality in brain dynamics and cardiovascular control’.

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