The Neurophysiological Biomarker Toolbox (NBT)

# Part IV: Granger-causality indices

In this part of the tutorial we will generate surrogate EEG data using the VAR model that we just fitted to our real EEG measurements. However, remember that scalp EEG data is affected by volume conduction effects. Consequently, a GC analysis of our real EEG is likely to result in a fully connected graph, i.e. in the absence of any relevant or interesting direction of information flow. But let's start checking whether this statement is really true.

### Two GC indices

In this tutorial we will use two popular linear indices that are based on the concept of Granger Causality:

These two indices are similar in the sense that both are expressed in the frequency domain and that both are _normalized_, i.e. they can take values from 0 to 1 with 0 meaning no flow and 1 meaning a maximal information flow.

### PDC/DTF Analysis of the scalp EEG data

A DTF and PDC analysis always start by fitting a VAR model to your data:

arfitObj = learn(var.arfit, eeg);

Then you can compute the DTF with the command:

dtfObj = compute(var.dtf, arfitObj, linspace(0,0.1,50));

Or the PDC with the command:

pdcObj = compute(var.pdc, arfitObj, linspace(0,0.1,50));

The last input argument to method compute() is an array of normalized frequencies at which the PDC and DTF will be computed. Remember that a normalized frequency is just the result of dividing a real frequency by the sampling rate. This means that the commands above will compute the PDC and DTF in a range of frequencies starting at 0 Hz and ending at $0.1\times (1000)=100$ Hz. Within this range the PDC and DTF values will be computed at 49 equidistant bins (that is what function linspace does).

### Plotting the PDC/DTF results

You can plot the result of the DTF and PDC analyses using these commands:

figure;
plot(dtfObj);
set(gcf, 'Name', 'DTF analysis of the scalp EEG data');
figure;
plot(pdcObj);
set(gcf, 'Name', 'PDC analysis of the scalp EEG data');

The interpretation of the figures is as we explained during the lecture. Cell $(i,j)$ of the figure corresponds to the flow in the direction $j\rightarrow i$. The vertical axis of each cell is PDC or DTF value (ranging from 0 to 1). The horizontal axis of each cell is normalized frequency. You should be able to zoom into a cell element by clicking on it. You can then return to the original figure by clicking on the figure again.

### Significance thresholds

Since both the PDC and DTF are just estimates, they are affected by random estimation errors. Therefore, it is necessary to assess how large a PDC or DTF value needs to in order to consider it significantly larger than 0. This can be done using method compute_significance():

pdcObj = compute_significance(pdcObj, surrogator.var, 0.05)
dtfObj = compute_significance(dtfObj, surrogator.var, 0.05)

where the last argument is the significance threshold. Computing the significance thresholds may take a few minutes. Once done you can display the thresholds by plotting again objects pdcObj and dtfObj:

figure;
plot(dtfObj);
set(gcf, 'Name', 'DTF analysis of the scalp EEG data');
figure;
plot(pdcObj);
set(gcf, 'Name', 'PDC analysis of the scalp EEG data');

which will generate the following two figures:

There is no point in trying to interpret the results of this analysis. The EEG data that you have been provided is just a very small chunk of data from a single subject and from just 7 (randomly chosen) EEG sensors (out of 257 sensors). Performing and interpreting the results of a real study is far beyond the scope of this tutorial. However, what you can conclude from the figures above is that there are significant flows at one frequency or another in almost any direction, i.e. the PDC (and even less the DTF) do not seem to be very specific. The probable reason for this is that information flows between scalp EEG signals is greatly affected by volume conduction effects.

### The difference between the PDC and the DTF?

The DTF and PDC analysis of the scalp EEG data revealed that the two indices produce different results. From that analysis is not easy to tell what the difference is. However, you can use the tools described in Part II of the tutorial to simulate VAR observations of known topology (e.g. a path-like topology). Then you can analyze this artificial observations in the same way as we did above with the scalp EEG data and draw some conclusions on the results provided by the PDC and the DTF.

##### Hint:

One of the two indices is mostly sensitive to direct flows of information while the other is also very sensitive to indirect flows. Which one is which? In a network with many connections between nodes, do you think that it is better to have an index that is sensitive to indirect flows, or is it better to be sensitive only to direct flows?

### Create an EEG-inspired VAR model

We will now create a VAR model that emulates the model that we fitted to the scalp EEG in the sense that the generated observations will look like EEG signals. However, in our artificial model we will modify the pattern of information flow between EEG signals so that it matches our wishes. This can be done with the commands:

varObj = dynamics.var('VarCoeffs', arfitObj.Coeffs);
varObj = set_topology(varObj, 'myTopo', 'randomize', false);

where 'myTopo' should be replaced by the topology that you want your artificial EEG data to have. You can use any of the topologies that we discussed in Part II of the tutorial, with the exception of the 'random' topology.

Now you should draw in a paper the information flow diagram underlying your synthetic EEG model. You can do this by inspecting the model coefficients of your model (varObj.VarCoeffs), as we described in Part II of the tutorial.

### Generate artificial EEG

Now you should generate some observations of the EEG-inspired VAR model that you just created:

[artificialEEG varObj] = generate(varObj, 10000);

Plot this artificial EEG and confirm that, contrary to the toy model that we used in Part II, these observations do look like real EEG:

eegplot(cell2mat(artificialEEG), 'srate', 1000);

In my case, the model observations look like this:

Now you can save your artificial EEG in a file called myeeg.mat:

artificialEEG = cell2mat(artificialEEG);
save myeeg.mat artificialEEG;

And finally, you can send this file to one of your fellow students and challenge him or her to guess the pattern of information flow underlying your artificial EEG data. For that, he or she will have to perform an analysis of your artificial EEG data in the same way as we did with the real EEG data.

##### Hint:

You are allowed to read the help of the constructor of dynamics.var objects:

help dynamics.var;

and build a more complicated VAR model (e.g. with noise, with non-equal innovation variances, etc…).

##### Hint:

You can send the file myeeg.mat either via e-mail, with a USB memory stick (ask your tutor for one) or via http://www.wetransfer.com.

### References

[1] Kaminski and Blinowska, 1991, A new method of the description of the information flow in the brain structures, Biol. Cybern. 65, 203-210, doi: 10.1007/BF00198091

[2] Baccala and Sameshima, 2001, Partial Directed Coherence: a new concept in neural structure determination, Biol. Cybern. 84, 463-474, doi: 10.1007/PL00007990