Recurrence Plots and Recurrence Quantification of EEG

These are the lecuture notes to a crashcourse introduction to the analysis of EEG data by means of recurrence quantification. The course was held in Vietri sul Mare (IT), May 2-4, 2012 within the timely (http://www.timely-cost.eu/) framework and supported by European Union.

Contents

Requirements

In order to work through the examples you will need Matlab and the CRPtoolbox which can be downloaded here: http://tocsy.pik-potsdam.de/CRPtoolbox/.

The toolbox is written and maintained by Norbert Marwan of the Potsdam Institute for Climate Impact Research.

Comments/Questions:

Are welcome and should be sent to: schinkel@physik.hu-berlin.de

Definition Recurrence plot

A recurrence plot (RP) is an advanced technique of nonlinear data analysis. It is a visualisation (or a graph) of a square matrix, in which the matrix elements correspond to those times at which a state of a dynamical system recurs (columns and rows correspond then to a certain pair of times). Techniqually, the RP reveals all the times when the phase space trajectory of the dynamical system visits roughly the same area in the phase space.

System Dynamics, Phase spaces and their reconstruction

To get an idea about phase spaces, we will start by investigating a prototypical system, the Lorenz System. We will first look at its phase space portrait and the reconstruct the phase from only a partial observation. The reconstruction of a systems phase space is often a requirement for the construction of recurrence plots from measured data.

Lorenz system

The Lorenz system in a common toy system in nonlinear dynamics. It is defined by the following equations:

$$\dot{x}~=~\sigma (y-x)$$

$$\dot{y}~=~\rho x-y-xz $$

$$\dot{z}~=~-\beta z + xy $$

with the ] parameters $\rho=28$, $$\sigma = 10$, $\beta = 8/3$ the system evolves around a strange attractor. The phase space portrait is the famous butterfly-like shape.

This function can easily be solved in Matlab. Save the function definition below to a file called lorenzEq.m and the use a built-in solver (ode45) to integrate the equations.

% function dx = lorenzEq(t,x)
%
% dx = zeros(3,1);
% dx(1) = 10 * (x(2) - x(1));
% dx(2) = 28 * x(1) - x(2) - x(1) * x(3);
% dx(3) = x(1) * x(2) - 8/3 * x(3);

Solving the Lorenz System

solve the equation

[t x] = ode45('lorenzEq',[0:.01:50],[5 5 5]);


% To get a head start we can use precomputed data
% stored on my website
url= 'http://dreeg.org/timely/lorenz.dat';
x = urlread(url);
% since it was stored as text, we have to make it to numbers
% and reshape it a little
x = str2num(x);
t = x(:,1); x = x(:,2:4);

% plot lorenz
plot(t,x);
xlabel('Time'), ylabel('x, y, z'),
title('Lorenz System')

% plot ps portrait
figure
plot3(x(:,1), x(:,2),x(:,3));
view(-173,42);
xlabel('x'), ylabel('y'), zlabel('z');

Phase space reconstruction

Now we reconstruct the system from only on component. This is the delay-embedding introduced by Takens

$$\vec{x}(i)~=~ (u_i,u_{i+\tau},\dots,u_{i+\tau(m-1)})$$

figure
tau = 22;
plot3(x(1:end-2*tau,1),x(1+tau:end-tau,1),x(1+2*tau:end,1));
view(2);
xlabel('x'), ylabel('y'), zlabel('z')

% The system is recognisable, though a little distorted, by the main,
% feature, that it evolves around to fix points restored. This is even
% more obvious if we look at the temporal dynamics.

Temporal dynamics of the Lorenz System

The interesting thing about the Lorenz system is its evolution over time. Lets look at that

% set up a figure
figure
hold on
view(-173,42);
xlim([-20 20]);
ylim([-30 30]);
zlim([0 50]);

%loop over the data to see the temporal evolution of the system
for i = 1:5:3000
	plot3(x(i:i+10,1),x(i:i+10,2),x(i:i+10,3));
	drawnow
end

% now we do the same for the reconstructed data
clf
hold on
xlim([-20 20]);
ylim([-20 20]);
view(2);
for i = 1:5:3000,
	plot3(x(i:i+10,1),x(i+tau:i+tau+10,1),x(i+2*tau:i+2*tau+10,3));
	drawnow
end

Parameter estimation

The required parameters can be estimated using functions in the CRPtoolbox

% The mi function computes the time-lagged mutual information function
% we will use the first minimum as an estimate for the embedding delay.
% It uses histograms, we use 10 bins for those, and compute for a maximal
% lag of 30.

nBins = 10;
maxDelay = 30;
mi(x(:,1),nBins,maxDelay);

% The FNN function implement a false nearest neighbours approach
% to estimate the embedding dimension, we go up to dimension 10
% and use the delay estimated above(22)
maxDim = 10;
tau = 22;
fnn(x(:,1),maxDim,tau);

Recurrence Plots

Lets now get to recurrence plot. We already have a chaotic system, but for comparison it would be nice to have some random and periodic system as well.

% construct noise and sine
sine = sin( 0:pi/100 : 5*pi );
noise = rand(1,length(sine));

"Standard" RPs

Now we compute the RPs for all three systems using the euclidean norm (numeric distance in phase space)

rpSine = crp(sine,2,1,.1,'eucl','norm');
rpNoise = crp(noise,2,1,.1,'eucl','norm');
rpLorenz = crp(x(:,1),3,22,.3,'eucl','norm');

use method: Euclidean norm
normalize data
use method: Euclidean norm
normalize data
use method: Euclidean norm
normalize data

Visual analysis

Visual inspection already reveals that typical system have a unique recurrence structure. In noise there is hardly any structure, periodical system have long continuous lines and chaotic systems are somewhere in between.

figure;imagesc(rpNoise);title('noise')
set(gca,'Ydir','normal')
colormap(flipud(gray(2)));xlabel('time'),ylabel('time');

figure;imagesc(rpSine);title('sine')
set(gca,'Ydir','normal')
colormap(flipud(gray(2)));xlabel('time'),ylabel('time');

figure;imagesc(rpLorenz);title('lorenz')
set(gca,'Ydir','normal')
colormap(flipud(gray(2)));xlabel('time'),ylabel('time');
xlim([1 1000])
ylim([1 1000])

Order Patterns Recurrence Plot

Instead of a numerical similarity, we can also use order patterns to look for recurrences in a symbolic space

oprpSine = crp(sine,2,1,.1,'op','norm');
oprpNoise = crp(noise,2,1,.1,'op','norm');
oprpLorenz = crp(x(:,1),3,22,.3,'op','norm');

figure;imagesc(oprpNoise);title('noise')
set(gca,'Ydir','normal')
colormap(flipud(gray(2)));xlabel('time'),ylabel('time');

figure;imagesc(oprpSine);title('sine')
set(gca,'Ydir','normal')
colormap(flipud(gray(2)));xlabel('time'),ylabel('time');

figure;imagesc(oprpLorenz);title('lorenz')
set(gca,'Ydir','normal')
colormap(flipud(gray(2)));xlabel('time'),ylabel('time');

use method: order pattern
normalize data
use method: order pattern
normalize data
use method: order pattern
normalize data

Patterns vs. Numbers

Again looking at the figures, what do you see? What can you learn for the analysis based on order patterns?

Recurrence Quantification

Recurrence Quantification can be used to further investigate the structures in an RP, which provides insights on the time dependent behaviour of the system. We will practice that with EEG data from a standard oddball-experiment (Sutton 1965)

EEG analysis using RQA

EEG/ERP analysis is used to investigate the brains response to given stimulus. The prime advantage it its high temporal resolution. In order to study the temporal dynamics, we will use the RQA in a windowed fashion ie. we extract small windows from the RP and quantify those.

% The data for a single subject P3 can be downloaded from my website.
% For convenience we will refer to the condition with P3 as "a" and that
% without as conditions "b"

% clean up workspace
clear
close all

% load data
urlA = 'http://dreeg.org/timely/c10.dat';
urlB = 'http://dreeg.org/timely/c11.dat';

a = str2num(urlread(urlA));
b = str2num(urlread(urlB));

% both matrices are 40x276, where the 1st dimension are the electrodes,
% the 2nd is time. The data is already averaged over trials.
% The EEG is sampled from -300ms pre-stimulus to 800 ms post-stimulus
% at 250Hz (4ms steps).

timeScale = -300:4:800;

% First we look at Voltage measurements at Pz (channel 27)
% which reveals a well-pronounced P300
plot(timeScale,a(27,:));
hold on;
plot(timeScale,b(27,:),'r');
legend('A','B');
xlabel('time in ms');
ylabel('Voltage');

Recurrence Plots of EEG data

Now we look at the RPs of the data. We will use the recurrence criterion introduced by Eckmann et al, fixed amount of nearest neighbours ('fan'). 

% Embedding dimension and delay 
dim = 3;
tau = 4;

% the threshold we use is .1 (10% nearest neighbours)
thresh = .1;

% construct RPs
rpA = crp(a(27,:),dim,tau,thresh,'fan');
rpB = crp(b(27,:),dim,tau,thresh,'fan');

% adjust the timeScale to compensate for embedding
tsNew = timeScale(1:size(rpA,1));

%and plot the RPs
figure;imagesc(rpA);
title('Condition A')
set(gca,'Ydir','normal')
set(gca,'XTickl',tsNew(get(gca,'XTick')))
colormap(flipud(gray(2)));
xlabel('time');
ylabel('time');

figure;
imagesc(rpB);
title('Condition B')
set(gca,'Ydir','normal')
set(gca,'XTickl',tsNew(get(gca,'XTick')))
colormap(flipud(gray(2)));
xlabel('time');
ylabel('time');



use method: fixed amount of nearest neighbours
normalize data
use method: fixed amount of nearest neighbours
normalize data

Look at the plots and inspect them for differences. Try different norms and thresholds (this is easiest done using the GUI 

crp(a(27,:);
crp(b(27,:);

Recurrence Quantification of EEG data

The plots already "look different". In order to quantify the differences, we compute the recurrence based complexity measures (ie RQA) over time in windows of 240ms with a shifting of 20ms.

win = 60;
step = 5;
% Try different norms and thresholds.
% We do this for all electrodes this time.
% Computation will take a little (about 2min
% or so). To not get too bored, we ask
% Matlab to tell us how long it takes
% for each electrode.

for iElec = 1:40
    tic,
	rqaA(iElec,:,:) = crqa(a(iElec,:),dim,tau,thresh,win,step,'fan','sil');
	rqaB(iElec,:,:) = crqa(b(iElec,:),dim,tau,thresh,win,step,'fan','sil');
    disp(sprintf('Electrode %d of %d took %2.2f seconds',iElec,40,toc));
end %iElec

Electrode 1 of 40 took 6.15 seconds
Electrode 2 of 40 took 5.99 seconds
Electrode 3 of 40 took 6.01 seconds
Electrode 4 of 40 took 6.00 seconds
Electrode 5 of 40 took 5.98 seconds
Electrode 6 of 40 took 6.00 seconds
Electrode 7 of 40 took 6.02 seconds
Electrode 8 of 40 took 6.01 seconds
Electrode 9 of 40 took 6.00 seconds
Electrode 10 of 40 took 5.99 seconds
Electrode 11 of 40 took 6.09 seconds
Electrode 12 of 40 took 6.05 seconds
Electrode 13 of 40 took 5.97 seconds
Electrode 14 of 40 took 6.00 seconds
Electrode 15 of 40 took 5.99 seconds
Electrode 16 of 40 took 6.00 seconds
Electrode 17 of 40 took 5.96 seconds
Electrode 18 of 40 took 5.99 seconds
Electrode 19 of 40 took 6.04 seconds
Electrode 20 of 40 took 6.01 seconds
Electrode 21 of 40 took 6.01 seconds
Electrode 22 of 40 took 6.00 seconds
Electrode 23 of 40 took 5.99 seconds
Electrode 24 of 40 took 6.03 seconds
Electrode 25 of 40 took 5.99 seconds
Electrode 26 of 40 took 6.03 seconds
Electrode 27 of 40 took 6.06 seconds
Electrode 28 of 40 took 6.02 seconds
Electrode 29 of 40 took 6.03 seconds
Electrode 30 of 40 took 6.02 seconds
Electrode 31 of 40 took 5.98 seconds
Electrode 32 of 40 took 6.00 seconds
Electrode 33 of 40 took 5.99 seconds
Electrode 34 of 40 took 6.00 seconds
Electrode 35 of 40 took 6.03 seconds
Electrode 36 of 40 took 5.98 seconds
Electrode 37 of 40 took 5.99 seconds
Electrode 38 of 40 took 5.98 seconds
Electrode 39 of 40 took 6.02 seconds
Electrode 40 of 40 took 6.01 seconds

Comparison of Conditions

rqa = crqa() returns DET in the second column (ie rqa(:,2) ), which is now shifted by one dimension since we did it for all electrodes. The measure DET is a measure for determinism (in the sense of predictability) in the the signal. Lets plot that for all electrodes in a surface plot.

figure;
surf(timeScale(1:step:size(rqaA,2)), 1:40,rqaA(:,1:step:end,2));
view(2);
colorbar;
axis tight
title('DET in condition A');
xlabel('time in ms');
ylabel('electrodes');

figure;
surf(timeScale(1:step:size(rqaB,2)), 1:40,rqaB(:,1:step:end,2));
view(2);
colorbar;
axis tight
title('DET in condition B');
xlabel('time in ms');
ylabel('electrodes');

This already reveals a higher level in DET in condition B, the condition that features a P300. To highlight that difference we now plot the difference

% compute difference
deltaDET = rqaB(:,:,2) - rqaA(:,:,2);

% and plot it
figure
surf(timeScale(1:step:size(deltaDET,2)), 1:40,deltaDET(:,1:step:end));
view(2);
colorbar;
axis tight
title('DET difference condition B-A');
xlabel('time in ms');
ylabel('electrodes');

Further reading/getting help

Plenty of information can be found on the net. A good starting point is http://www.recurrence-plot.tk, which has plenty of information about the method as well as a comprehensive bibliography and a forum.

Exercises:

1) Analyse the logistic map $\qquad x_{n+1} = r x_n (1-x_n)$ Begin with an interval for r from [3.5 4]. Then proceed by studying $r$ in [3 4]. Can you detect transitions ? Can you relate those to the bifurcation diagram of the logitic map

% function out = logistic(x0,r,n)
%   out(1) = x0;
%   for i = 1:n-1
%       out(i+1) = r*out(i)*(1-out(i));
%   end

2) Further study the EEG data provided. Use the software on my website http://dreeg.org to estimate the confidence bounds of RQA measures. (cf. Schinkel et al, 2009)

3) Try to construct networks based on joint recurrences (cf. Schinkel et al, 2011).


Published with MATLAB® 7.11