A SIMPLE TUTORIAL TO START WITH FUNCTIONAL NETWORKS.
Here we are! The main idea is to know how we can build up a functional brain network from a EEG recording database, then visualize some network features in the brain network. So let's star!
Contents
- Author
- Brief Data Description
- An EEG and its power.
- Band decompotition of one signal
- Gahtering dataset of 20 subject in condition 2 (mild sedation)
- Filtering alpha band the previous dataset
- Covariance
- Pearson
- Phase Locking Value (PLV)
- Magnitude-Squared Coherence
- Ordinal Synchronization (A new information-based correlation)
- Fully and not fully connected graph
- Some Networks Features
- Rescaling network features for plotting networks
- Plotting a Functional Network
- Acknowledgemetns
Author
Johann H. Martínez
Universidad de los Andes. Depto. Biomedical Engineering (2019) Bogotá. Colombia
https://johemart.wixsite.com/neurocomplexity
Brief Data Description
20 healthy individuals who were given the sedative propofol.
clearvars; clc; close all; % command for clear screen, command windows and clossing all open figures load ALLEEG.mat; % % 91 channels , 2500 timpepoints , 20 subject x 4 conditions = 80 essays load('-mat','02-2010-anest 20100210 135.003.set') % % extra info of dataset, e.g., time in ms, sampling rate in Hz and bands.
An EEG and its power.
Here, it is a single time series, then we plot a power of 91 signals of the same subject to see diferent power peaks. Note how the second subplot highlight the frequency bands.
ch = 1; % Taking channel one cd = 2; % taking second condition x = data(ch,:,cd); % this is the resulting time series figure; subplot(121) plot(EEG.times*1e-3, x) grid on xlabel('Time (s)') ylabel('Voltage (mV)') title('A signal in time domain') subplot(122) plot(EEG.freqs, EEG.spectra(:, :, 1)) xlim([0 45]) xlabel('frequency (Hz)') ylabel('Power (dB)') title('Powers of 91 channels-one subject') ylim([-40 30]) grid on hold plot([EEG.freqwin(2) EEG.freqwin(2)],[-40 30], '--k') plot([EEG.freqwin(3) EEG.freqwin(3)],[-40 30], '--k') plot([EEG.freqwin(4) EEG.freqwin(4)],[-40 30], '--k') plot([EEG.freqwin(5) EEG.freqwin(5)],[-40 30], '--k') plot([EEG.freqwin(6) EEG.freqwin(6)],[-40 30], '--k') text(2,27.5, '\delta', 'FontSize', 20) text(6,27.5, '\theta', 'FontSize', 20) text(10,27.5, '\alpha', 'FontSize', 20) text(17,27.5, '\beta', 'FontSize', 20) text(35,27.5, '\gamma', 'FontSize', 20) set(gcf, 'units','normalized','outerposition',[0 0 1 1]) % EXPANDING FIGURE ON SCREEN
Current plot held
Band decompotition of one signal
The signal is now filtered under different frequency bands. Here we can observe the broad band signal and its recontructions from 5 rhythms. The function " bandpass2.m " takes a time series " X ", the sampling rate, here as " EEG.srate ", and the limits of the frequency band.
xb1 = bandpass2(0, 4, EEG.srate, x); % signal filtered in delta band xb2 = bandpass2(4, 8, EEG.srate, x); % signal filtered in theta band xb3 = bandpass2(8, 15, EEG.srate, x); % signal filtered in alfa band xb4 = bandpass2(15, 25, EEG.srate, x); % signal filtered in beta band xb5 = bandpass2(25, 45, EEG.srate, x); % signal filtered in gamma band figure; subplot(611) plot(EEG.times*1e-3, x) ylabel('BB') title('Voltage band decomposition') subplot(612) plot(EEG.times*1e-3, xb1) ylabel('V-\delta') subplot(613) plot(EEG.times*1e-3, xb2) ylabel('V-\theta') subplot(614) plot(EEG.times*1e-3, xb3) ylabel('V-\alpha') subplot(615) plot(EEG.times*1e-3, xb4) ylabel('V-\beta') subplot(616) plot(EEG.times*1e-3, xb5) xlabel('Time (s)') ylabel('V-\gamma') set(gcf, 'units','normalized','outerposition',[0 0 1 1]) %EXPANDING FIGURE ON SCREEN
Gahtering dataset of 20 subject in condition 2 (mild sedation)
We state " cd = 2 " as the variable to choose the second condition. Observing the file " datainfo ", we detect that second column declares the level of sedation (condition) and the data organization (same condition each 4 files). Dataset fieds are: filename | sedation level (1=baseline, 2=mild, 3=moderate, 4=recovery) % propofol concentration in plasma (microgram/litre) | correct responses in a task
load datainfo.mat; % file to reveal the order at which date has been saved cd = 2; subper_cd = cd:4:80; % by observing "datainfo.mat" we know the same condition is saved each four files subs_cd2 = data(:,:,subper_cd); % 91 time series of 2500 points for 20 subjects under mild sedation
Filtering alpha band the previous dataset
As we observed the alpha band as the one of higher power, we filter previous dataset of mild sedation in alpha frequency band. Variable " subs_cd2_al " contains the 91 times series for 20 subjects in mild sedation for alpha.
subs_cd2_al = zeros(size(subs_cd2)); for su = 1 : 20 for ch = 1 : 91 subs_cd2_al(ch,:,su) = bandpass2(8, 15, EEG.srate, subs_cd2(ch, :, su)); end end
Covariance
Covariance gives the level of variability of random variables respect their means.
cov_cd2 = cov(subs_cd2_al(:,:,1)'); cov_cd2 = abs(cov_cd2 - diag(diag(cov_cd2))); figure; imagesc(cov_cd2); axis square; cb = colorbar; cb.Label.String = 'Covariance level'; title('|Covariance|') ylabel('nodes'); xlabel('nodes')
Pearson
Linear correlaton between two random variables. In other words, how close the interdependency is near to a linear relation. It's also known to be the Covariance normalized as well as the of the Cross-Correlation function.
pea_cd2 = corrcoef(subs_cd2_al(:,:,1)'); pea_cd2 = abs(pea_cd2 - diag(diag(pea_cd2))); figure; imagesc(pea_cd2); axis square; colorbar; title('|Pearson|')
Phase Locking Value (PLV)
How well waves' phases are locked in a phase synchrony. The main idea is to take the phases of waveforms by taking into account the real part of the Hilbert Transform of signals. Then applying the equation.
Here we take phases for subject one. Variable " phis " is then the phases of 91 signals of subject one.
sujeto = 1; phis = angle(hilbert(subs_cd2_al(:,:,sujeto)'));
Now computing PLV over phases "phis" of subject one.
pairs = nchoosek(1:91,2); % all combinations of 91 channels in groups of two plv = zeros(91,91); for pr = 1 : length(pairs) plv(pairs(pr,1), pairs(pr,2)) = abs(sum(exp(1i* ( phis(:,pairs(pr,1)) - phis(:,pairs(pr,2)) ))))/length(phis); end plv = plv + plv'; % we take its transpose to create the connectivity matrix figure imagesc(plv);colorbar; axis square; colorbar title('PLV')
Magnitude-Squared Coherence
Linear correlation between power espectrum in a frequency f . It is sensitive to both phase and amplitude. We take the average over the segments of the FFT.
We take two signals of subject one and visualize its coherence. When computing coherence, e.g., in alpha band, we will average the coherence points in the limits of such band (Here, they are red dots), and that mean value is called as the coherence for two signals in alpha band.
sujeto = 1; [coh, f] = mscohere(subs_cd2(1, :, sujeto), subs_cd2(2 ,:, sujeto), hanning(200), 25, 200, EEG.srate); figure plot(f,coh, '-') hold plot(f(f>=EEG.freqwin(3) & f<=EEG.freqwin(4)), coh(f>=EEG.freqwin(3) & f<=EEG.freqwin(4)), 'ro') plot([EEG.freqwin(3) EEG.freqwin(3)], [min(coh) 1], '--k') plot([EEG.freqwin(4) EEG.freqwin(4)], [min(coh) 1], '--k') legend('Coherence', '\alpha Coherence') xlabel('Frequency (Hz)') ylabel('Coherence') text(10,0.98, '\alpha', 'FontSize', 20) set(gcf, 'units','normalized','outerposition',[0 0 1 1]) %EXPANDING FIGURE ON SCREEN
Current plot held
Now we will take all signals of subject one in order to compute the coherence of all pair of signals in alpha band.
% % % sujeto = 1; % % % pairs = nchoosek(1:91,2); % % % cxy = zeros(91,91); % % % for pr = 1 : length(pairs) % % % [coh, f] = mscohere(subs_cd2(pairs(pr,1), :, sujeto), subs_cd2(pairs(pr,2) ,:, sujeto), hanning(200), 25, 200, EEG.srate); % % % cxy(pairs(pr,1), pairs(pr,2)) = mean(coh(f>=EEG.freqwin(3) & f<=EEG.freqwin(4))); % % % end % % % cxy = cxy + cxy'; % % % % % % figure % % % imagesc(cxy); % % % axis square; colorbar; title('Coherence')
Ordinal Synchronization (A new information-based correlation)
Fast, robust-to noise tool to assess synchronization at higher orders D. It was recently published in Physica A https://www.sciencedirect.com/science/article/abs/pii/S0960077918309081. When computing it, we take the average over M segments of ordinal vectors of dimension D.
% % % sujeto = 1; % % % pairs = nchoosek(1:91,2); % % % osn = zeros(91,91); % % % for pr = 1 : length(pairs) % % % osn(pairs(pr,1), pairs(pr,2)) = f_op_syn(subs_cd2_al(pairs(pr,1), :, sujeto), subs_cd2_al(pairs(pr,2), :, sujeto), 7, 'cons'); % % % end % % % osn = abs(osn + osn'); % % % % % % figure; % % % imagesc(osn); % % % axis square; colorbar % % % title('Ordinal Synchronization')
Fully and not fully connected graph
By taking a fully connected graph, one takes into account spurious correlations and we don't want that, so we need to threshold our matrix. Bad news. There is not an exact threshold. Good news! There are many ways to thresholding up a matrix.
Here we keep the of the weighted links of the Pearson correation matrix, always keeping the network conectedness. Here we make use of functions " threshold_proportional.m " and " isconnected " to threshold the matrix and testing its connectedness, respectively.
W = threshold_proportional(pea_cd2, 0.2); % function for threshold a matrix
isconnected(W); % function for testing its connectedness
figure subplot(221) imagesc(pea_cd2); colorbar; axis square title('Fully connected matrix') subplot(222) plot(graph(pea_cd2)) axis off title('What an ugly network!') subplot(223) W = threshold_proportional(pea_cd2, 0.2); % look at this! "W" is now our matrix with threshold isconnected(W); % if this value is one, then "W" is connected imagesc(W); colorbar; axis square title('Thresholded matrix') subplot(224) plot(graph(W)) axis off title('A much better network to work...') set(gcf, 'units','normalized','outerposition',[0 0 1 1]) %EXPANDING FIGURE ON SCREEN
Some Networks Features
Network features are calculated on the matrix representation of the graph. Here let's consider the previous pearson correlation matrix "_W_" to play with.
Strength is one of the first centrality measures. It compute the " strenght " of a node by summing up its weights. Considering a conectivity matrix W , its strenght reads for:
st = sum(W,2); figure; hold stem(st) [st_max, node_st_max] = max(st); stem(node_st_max, st_max, 'r') ylabel('Strength'); xlabel('nodes')
Current plot held
We also may obtain another information about network by computing its Probability Density Function (PDF) of finding a specific strenght
[N,binlimit] = histcounts(st, linspace(min(st),max(st),50)); % N is the frecuency of nodes along several bin limits figure hold; bar(binlimit(1,2:end), N/numel(N)) ylabel('P(k_i=k)'); xlabel('k')
Current plot held
Clustering is another basic measure of node centrality. It takes into account the proportion of nodes's triangles formed for each node respect to the total number of triangles. In other words:
cw = clustering_coef_wu(W); figure plot(cw) hold; [cw_max, node_cw_max] = max(cw); stem(node_cw_max, cw_max, 'or') ylabel('Weighted Clustering'); xlabel('nodes')
Current plot held
Eigenvector Centrality is a centrality measure that give information about whether a node is well connected or not, understanding "well connected" by that node that is connected with hubs. Mathematically speaking, this vector centrality is the eigenvector associated to the first eigenvalue when solving:
ec = eigenvector_centrality_und(W); figure plot(ec) hold; [ec_max, node_ec_max] = max(ec); stem(node_ec_max, ec_max, 'or') ylabel('eigenvector centrality'); xlabel('nodes')
Current plot held
Average Shortest Path as its name indicates, it's the average shortest path to connect two nodes randomly choosen in the network. In weighted graphs: the weighted the link, the closer two node are. We take the inverse of the weights of matrix W to build up the connection-length matrix ( L ). We apply the Dijkstra's algorithm (the most accepted) to find al shortest paths among nodes in a matrix called ( D ). The average shortest path is then the mean value of this resulting matrix.
L = weight_conversion(W, 'lengths'); % taking the inverse of matrix elements D = distance_wei(L); % applying Dijkstra's algorithm sp = charpath(D); % averaging shortest paths' matrix figure plot(1,sp, '*k') hold; boxplot(D(D~=0)); legend(horzcat('sp = ', num2str(sp))) ylabel('Shortest Path')
Current plot held
Rescaling network features for plotting networks
First, you have to define what network feature you want highlight by plotting a functional network (a node parameter, a link feature?). Then you have to rescale that feature respect to numerical figure's parameters of the language you are using to. In this example we will plot the Strenght and link weights.
Although there are several ways for scaling a collection of numbers, for simplicity here we adopt the classical linear scaling:
Where a and b are the rescaling limits.
grafo = graph(W); % converting matrix to graph list of edges links = table2array(grafo.Edges); % converting the graph object into a table wij = links(:,3); % taking only link weigths a = 1; b = 3; % range for rescaling limits wij = ( (wij - min(wij)).*(b-a) )/(max(wij) - min(wij)); % rescaling formula wij = wij+0.5; rgb = squeeze(double2rgb(wij, colormap(jet))); close % mapping values to RGB colors with pallete JET a = 1; b = 3; % range for rescaling limits wij = ( (wij - min(wij)).*(b-a) )/(max(wij) - min(wij)); % rescaling formula wij = wij+0.5;
Plotting a Functional Network
The main idea is to superimpose different figures in a same canvas. Here we plot links, and then we superimpose nodes' plot.
figure hold on; for ll = 1 : length(links) %% plotting links line([EEG.chanlocs(links(ll,1)).X, EEG.chanlocs(links(ll,2)).X], [EEG.chanlocs(links(ll,1)).Y, EEG.chanlocs(links(ll,2)).Y], 'LineWidth', wij(ll), 'Color', rgb(ll,:) ); end colormap(jet); % colorbar for nn = 1 : 91 %% plotting nodes plot(EEG.chanlocs(nn).X, EEG.chanlocs(nn).Y, 'ro', 'MarkerSize', st(nn), 'MarkerFaceColor', 'y') text(EEG.chanlocs(nn).X+0.5, EEG.chanlocs(nn).Y, EEG.chanlocs(nn).labels, 'FontSize', 11); end view([270 90]); %% rotating to vertical view axis equal; axis off; set(gcf, 'units','normalized','outerposition',[0 0 1 1]) %EXPANDING FIGURE ON SCREEN
Acknowledgemetns
I would like to thank to my students and colleagues for comments in the preparation of this tutorial. We used some functions from the following tookits:
Brin Connectivity Toolbox by Rubinov M and Sporns O. https://sites.google.com/site/bctnet/
Matlab Tools for Network Analysis by MIT http://strategic.mit.edu/downloads.php?page=matlab_networks
function "double2rgb.m" by David Legland https://github.com/mattools/matImage/blob/master/matImage/imFilters/double2rgb.m
Dataset were extracted from the paper by Chennu et all PlosOne (2016) https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004669