Data preprocessing and most of the following steps of the data analysis were done using the fieldtrip toolbox (F. C. Donders Centre for Cognitive Neuroimaging: http://www.ru.nl/fcdonders/fieldtrip). First, all data sets were down-sampled to 600 Hz and cut into epochs of two seconds and those epochs containing blinks or muscle artifacts were excluded from further analysis based on visual inspection. Second, an independent component analysis (ICA) was calculated for each individual data set to identify components that reflect the heart-beat and these components were rejected from the data (using the logisitic infomax ICA algorithm implemented in eeglab: http://sccn.ucsd.edu/eeglab/). After artifact correction, 90 trials (i.e. 180 seconds in total) were selected randomly from the remaining trials and used for the following analyses. This selection was done to keep the number of trials constant across all subjects. The number of 90 trials reflects a trade-off between cleaning the data from noisy events as much as possible and still having enough data to calculate the autoregressive model.

#### Step 2: Partial directed coherence

For each subject, we computed partial-directed coherence (PDC) for the full set of voxels [

19,

20]. Partial- directed coherence is a measure of effective coupling that captures the direction of the information-transfer between the given voxels. Thus, with a set of N voxels, we get a total of NxN PDC-values for each subject that reflects for each pair of voxels the effective coupling in both directions. This approach is based on multivariate autoregressive (MVAR) modeling that integrates temporal and spatial information. Here, we model for each voxel the influence by all other voxels for a given time-range. The model order

*p* defines this time range of the autoregressive process and describes how many time points - back in time - are used for the modeling the current value. In the univariate case this can be written as

whereby

*y(t)* denotes the predicted value at time-point

*t*,

*a(1), a(2),...a(p)* determine the regression coefficient and

*x(t)* is called the

*innovation process* which equals the difference between the actual value at time

*t* and the estimation of

*y(t)* based on the linear combination of the previous time points

*y(t-1), y(t-2),... y(t-p)* [

43]. In order to find the optimal model parameter

*P* we calculated the Schwarz Bayesian Criterion (SBC) [

44] for model orders from 2 - 20. On average over the whole sample, the minimum of the SBC function was located at

*P* = 6 which was then taken as the model order for all subjects. For estimation of the autoregressive parameters we used the Vieira-Morf algorithm [

45] implemented in the biosig toolbox (

http://www.biosig.sf.net, version 2.12) which has been found to provide the most accurate estimates [

43]. The matrix of autoregressive coefficients in the multivariate case can be written as

where the coefficients *aij* represent the linear interaction between voxel *i* onto voxel *j* for a given time lag *k*.

Partial Directed Coherence is a statistical measure that is related to the concept of Granger Causality [

46] and is able to detect asymmetric coupling between the compared voxels for a given frequency range. Here we investigated the frequency range of 2 to 100 Hz (steps of 2 Hz). In order to reveal the spectral properties, the autoregressive coefficients are transformed into the frequency domain by

with

representing the matrix of the frequency-transformed autoregressive coefficents,

*I* being the identity matrix and

*f*
_{
s
}being the sampling frequency.

With denoting the

*i, j*-th element of the relative coupling strength from voxel

*j* to voxel

*i* at a given frequency

*f*, the directed information flow from

*j* to

*i* can be written by

The superscript *H* denotes the Hermetian transpose which is found by taking the complex conjugate of each entry of the standard matrix transpose. Thus, the PDC value π ii *(f)* indicates how much the activity of voxel *i* depends on its own past at a given frequency. The value π ij *(f)* denotes how much the frequency-specific activity of voxel *j* depends on voxel *i*. The PDC estimators were calculated using functions implemented in the biosig toolbox (http://www.biosig.sf.net, version 2.12).

To the best of our knowledge, there is no established way of calculating the statistical significance of the PDC estimators. Thus, we used a permutation approach to estimate thresholds for significant coupling between pairs of voxels (couplings of one voxel with itself were excluded from the analysis). Therefore, the following steps 1) to 3) were repeated 1,000 times:

1) Shuffle the matrix A of the autoregressive coefficients pseudo-randomly. This was done the following way: The matrix A is a square matrix with 326 rows and 326 columns. Firstly, we generated a vector with random numbers between 1 and 326. Secondly, the columns were shuffled according to the random vector. Thirdly, the rows were shuffled according the same random vector. This shuffling procedure was repeated for all model orders.

2) Calculate the PDC estimators in the way that was described above.

3) Determine the 99%-percentile of the PDC estimator for each frequency and save it. The 99%-percentile was used instead of the maximum to reduce the influence of the self-reflective coefficients (voxel *i* with itself) which are much higher and were not part of this analysis anyway.

The maxima over the 1,000 permutations was used as a threshold of significance for each frequency bin. Thresholds were calculated for each participant individually.