From 7580e72a4d05ed7815abfac572459966ad74c411 Mon Sep 17 00:00:00 2001 From: Yacine Mahdid Date: Tue, 21 Apr 2020 19:49:20 -0400 Subject: [PATCH] Add hub location to NeuroAlgo (#59) * add binary hub location function in the graph theory folder * document some more the binary_hub_location * add test for hub location * Fix up the already existing na_hub_location * add more helpful help message * add documentation on the readme.md about hub location --- MATLAB/example/example_analysis.m | 5 +- MATLAB/example/test_hub_location.m | 47 ++-------- .../graph_theory/binary_hub_location.m | 88 ++++++++----------- MATLAB/plot/make_video_hub_location.m | 25 ++---- MATLAB/src/na_hub_location.m | 58 +++++++++--- README.md | 15 ++++ 6 files changed, 121 insertions(+), 117 deletions(-) diff --git a/MATLAB/example/example_analysis.m b/MATLAB/example/example_analysis.m index bbd661ad..aadab537 100644 --- a/MATLAB/example/example_analysis.m +++ b/MATLAB/example/example_analysis.m @@ -27,6 +27,9 @@ recording class take a look at the /source folder %} + + + % Spectral Power window_size = 10; time_bandwith_product = 2; @@ -60,7 +63,7 @@ p_value = 0.05; % the p value to make our test on threshold = 0.10; % This is the threshold at which we binarize the graph step_size = 10; -%result_hl = na_hub_location(recording, frequency_band, window_size, step_size, number_surrogate, p_value, threshold); +result_hl = na_hub_location(recording, frequency_band, window_size, step_size, number_surrogate, p_value, threshold); % Permutation Entropy (PE) frequency_band = [7 13]; % This is in Hz diff --git a/MATLAB/example/test_hub_location.m b/MATLAB/example/test_hub_location.m index 23f91b41..31ec3a3b 100644 --- a/MATLAB/example/test_hub_location.m +++ b/MATLAB/example/test_hub_location.m @@ -16,57 +16,26 @@ step_size = 0.1; result_wpli = na_wpli(recording, frequency_band, window_size, step_size, number_surrogate, p_value); -figure; -subplot(2,1,1) -imagesc(squeeze(mean(result_wpli.data.wpli))); -colorbar -colormap jet; -title(strcat("Average Participant at ", string(number_surrogate), " surrogates")); -subplot(2,1,2); -imagesc(squeeze(result_wpli.data.wpli(5,:,:))); -colormap jet; -colorbar; -title(strcat("Single participant #",string(15)," at ", string(number_surrogate), " surrogates")); - -% Calculating Hub Location on the average wpli +% Calculating Hub Location on the each window of wpli channels_location = result_wpli.metadata.channels_location; % calculate hub location t_level_wpli = 0.2; -t_level_hub = 0.10; wpli = result_wpli.data.wpli; [num_window, num_channels,~] = size(wpli); -degree_map = zeros(num_window, num_channels); +% Here we are iterating on each window of the wpli matrices, binarize them +% and then running the binary_hub_location it +hub_map = zeros(num_window, num_channels); median_location = zeros(1,num_window); max_location = zeros(1,num_window); for i = 1:num_window % binarized the wpli matrix b_wpli = binarize_matrix(threshold_matrix(squeeze(wpli(i,:,:)), t_level_wpli)); - [normalized_location,previous_location, channels_degree] = binary_hub_location(b_wpli, channels_location, t_level_hub); - degree_map(i,:) = channels_degree; - median_location(i) = normalized_location; - max_location(i) = previous_location; + [hub_location, weights] = binary_hub_location(b_wpli, channels_location); + hub_map(i,:) = weights; end +% We are showing this stack of weights map into a small movie filename = 'hub_location_technique_comparison'; -make_video_hub_location(filename, degree_map, median_location, max_location, channels_location, step_size) - - -function plot_hub(normalized_location, previous_location, channels_degree, channels_location) - % Plotting the hub location - figure; - subplot(1,3,1) - topoplot(channels_degree,channels_location,'maplimits','absmax', 'electrodes', 'off'); - colormap('cool'); - colorbar; - title('Degrees') - subplot(1,3,2); - bar(normalized_location); - ylim([0 1]) - title('Current Normalized Location (median)'); - subplot(1,3,3); - bar(previous_location); - ylim([0 1]) - title('Previous Normalized Location (max)'); -end +make_video_hub_location(filename, hub_map, channels_location, step_size) diff --git a/MATLAB/features/graph_theory/binary_hub_location.m b/MATLAB/features/graph_theory/binary_hub_location.m index 270189f7..9eeff78d 100644 --- a/MATLAB/features/graph_theory/binary_hub_location.m +++ b/MATLAB/features/graph_theory/binary_hub_location.m @@ -1,71 +1,59 @@ -function [median_degree_location,max_degree_location, channels_degree] = binary_hub_location(b_matrix, channels_location, t_level) -%HUB_LOCATION choose from the the channels the channel with the highest -%degree -% Input: -% b_matrix: a binary undirected matrix -% t_level: what top percentage represent a hub in the brain (what is the -% definition of a hub) - - num_element = length(b_matrix); - - - %% Caculate the unweighted degree of the network - channels_degree = degrees_und(b_matrix); - - %% Calculate previous location - [~, channel_index] = max(channels_degree); - max_degree_location = threshold_anterior_posterior(channel_index, channels_location); +function [hub_location, weights] = binary_hub_location(b_wpli, location) +%BETWEENESS_HUB_LOCATION select a channel which is the highest hub based on +%betweeness centrality and degree +% input: +% b_wpli: binary matrix +% location: 3d channels location +% output: +% hub_location: This is a number between 0 and 1, where 0 is fully +% posterior and 1 is fully anterior +% weights: this is a an array containing weights of each of the channel in +% the order of the location structure - %% Threshold the degree to keep only t - sorted_degree = sort(channels_degree); - t_index = floor(num_element*(1 - t_level)) + 1; - t_element = sorted_degree(t_index); + %% 1.Calculate the degree for each electrode. + degrees = degrees_und(b_wpli); + norm_degree = (degrees - mean(degrees)) / std(degrees); + a_degree = 1.0; + %% 2. Calculate the betweeness centrality for each electrode. + bc = betweenness_bin(b_wpli); + norm_bc = (bc - mean(bc)) / std(bc); + a_bc = 1.0; - %% Threshold the degrees - t_channels_degree = channels_degree; - t_channels_degree(t_channels_degree < t_element) = 0; - % Get only the hub degree - non_zero_hub_degree = []; - for i = 1:num_element - current_degree = t_channels_degree(i); - if(current_degree > 0) - non_zero_hub_degree(end+1) = current_degree; - end - end - % get the median value of the non-zero hub degree - % to detect which region should be returned (we don't look at extreme - % values) - m_hub_degree = find_hub(non_zero_hub_degree); - disp("The degree of the selected hub:"); - disp(m_hub_degree) - m_hub_degree_index = find(channels_degree == m_hub_degree); - median_degree_location = threshold_anterior_posterior(m_hub_degree_index, channels_location); -end + %% 3. Combine the two Weightsmetric (here we assume equal weight on both the degree and the betweeness centrality) + weights = a_degree*norm_degree + a_bc*norm_bc; + + %% Obtain a metric for the channel that is most likely the hub epicenter + [~, channel_index] = max(weights); + hub_location = threshold_anterior_posterior(channel_index, location); -% This will try to find the middle degree value not looking at extreme -% values -function [middle_value] = find_hub(hubs) - sorted_hub_degree = sort(hubs); - middline_index = floor(length(sorted_hub_degree)/2); - middle_value = sorted_hub_degree(middline_index); end function [normalized_value] = threshold_anterior_posterior(index,channels_location) -%THRESHOLD_ANTERIOR_POSTERIOR Summary of this function goes here -% Detailed explanation goes here +%THRESHOLD_ANTERIOR_POSTERIOR This function will squash the channels in the +%direction of the anterior-posterior line and will set the most posterior +%index to 0 and the anterior most anterior index to 1 +% To do so it takes the index of the channel we want to normalize in the +% posterior-anterior direction and the channels location for where it +% came from. It then will find the value this channel should have. +% input: +% index: index of the channel to normalize between 0 and 1 +% channels_location: 3D channel data structure + % Get the X index of the current channels current_x = channels_location(index).X; + % Accumulate all the Xs index for the channels locations all_x = zeros(1,length(channels_location)); for i = 1:length(channels_location) all_x(i) = channels_location(i).X; end + % Normalize the current value of x between the min coordinate we have + % in the headset and the maximum min_x = min(all_x); max_x = max(all_x); - normalized_value = (current_x - min_x)/(max_x - min_x); end diff --git a/MATLAB/plot/make_video_hub_location.m b/MATLAB/plot/make_video_hub_location.m index 1a5f60b0..5124b6ce 100644 --- a/MATLAB/plot/make_video_hub_location.m +++ b/MATLAB/plot/make_video_hub_location.m @@ -1,8 +1,8 @@ -function make_video_hub_location(filename, degree_map, median_location, max_location, channels_location, step_size) +function make_video_hub_location(filename, hub_map, channels_location, step_size) %make_video_functional_connectivity Summary of this function goes here % Detailed explanation goes here - [number_frames, number_channels] = size(degree_map); - all_edges = reshape(degree_map, [1 (number_frames*number_channels)]); + [number_frames, number_channels] = size(hub_map); + all_edges = reshape(hub_map, [1 (number_frames*number_channels)]); min_color = min(all_edges); max_color = max(all_edges); @@ -11,31 +11,22 @@ function make_video_hub_location(filename, degree_map, median_location, max_loca open(video); %open the file for writing for i=1:number_frames %where N is the number of images disp(strcat("Making video: ", string(i/number_frames)," %")); - channels_degree = squeeze(degree_map(i,:)); + hub = squeeze(hub_map(i,:)); % Create the figure - fc_figure = plot_hub( channels_degree, median_location(i), max_location(i), channels_location, min_color, max_color); + fc_figure = plot_hub( hub, channels_location, min_color, max_color); writeVideo(video,getframe(fc_figure)); %write the image to file delete(fc_figure) end close(video); %close the file end -function [figure_handle] = plot_hub( channels_degree, median_location, max_location, channels_location, min_color, max_color) +function [figure_handle] = plot_hub( hub_map, channels_location, min_color, max_color) % Plotting the hub location figure_handle = figure('visible','off'); - subplot(1,3,1) - topoplot(channels_degree,channels_location,'maplimits','absmax', 'electrodes', 'off'); + topoplot(hub_map,channels_location,'maplimits','absmax', 'electrodes', 'off'); colormap('cool'); colorbar; caxis([min_color max_color]) - title('Degrees') - subplot(1,3,2); - bar(median_location); - ylim([0 1]) - title('Median Location'); - subplot(1,3,3); - bar(max_location); - ylim([0 1]) - title('Max Location'); + title('Hub Location') end \ No newline at end of file diff --git a/MATLAB/src/na_hub_location.m b/MATLAB/src/na_hub_location.m index 7ae90aa3..81770c90 100644 --- a/MATLAB/src/na_hub_location.m +++ b/MATLAB/src/na_hub_location.m @@ -1,6 +1,41 @@ function [result] = na_hub_location(recording, frequency_band, window_size, step_size, number_surrogate, p_value ,threshold) -%NA_HUB_LOCATION Summary of this function goes here -% Detailed explanation goes here +%NA_HUB_LOCATION will calculate hub locations +%(as defined by betweeness-centrality and degree) for the whole recording +%of EEG data and store it inside a result data structure. +% +% this function will take a EEG data bundled inside an recording instance +% and will calculate wPLI on segment of the data. It will then binarize +% these matrices of connectivity before calculating the binary hub +% location. It will output both the X-coordinate position of the hub +% along the anterior-posterior axis and the weights calculated to defined +% a hub. The hub here is defined as the maximum values in the map formed +% by doing degree+ betweness_centrality. +% +% input: +% recording: a EEG data bundled into a recording instance +% frequency_band: the frequency band at which to filter the data +% window_size: the window size onto which to cut the segment of data +% step_size: the step size at which to jump from one window to another if +% its the same as the window size it will be a jumpy windowing. +% number_surrogate: the number of surrogate wPLI matrices to make the +% null distribution for the corrected wPLI calculation +% p_value: p value at which to say that a connection is significant +% threshold: from 0 to 1 it's the amount of top connection we want to +% keep in the wPLI->binary_wpli. +% +% output: +% result: a datastructure containing all the parameters information, +% metadata as well as the output of the hub location. +% +% example usage: +% recording = load_set('name_of_data.set',path_to_data); +% frequency_band = [7 13]; +% window_size = 10; +% number_surrogate = 10; +% p_value = 0.05; +% threshold = 0.10; +% step_size = 10; +% result_hl = na_hub_location(recording, frequency_band, window_size, step_size, number_surrogate, p_value, threshold); %% Getting the configuration configuration = get_configuration(); @@ -21,25 +56,28 @@ % Here we init the sliding window slicing recording = recording.init_sliding_window(window_size, step_size); number_window = recording.max_number_window; + location = recording.channels_location; %% Calculation on the windowed segments + result.data.wpli = zeros(number_window, recording.number_channels, recording.number_channels); result.data.hub_index = zeros(1,number_window); - result.data.hub_degree = zeros(1,number_window); - result.data.hub_relative_position = zeros(1,number_window); - result.data.graph = zeros(number_window, recording.number_channels, recording.number_channels); + result.data.hub_weights = zeros(1,recording.number_channels); for i = 1:number_window print_message(strcat("Hub Location at window: ",string(i)," of ", string(number_window)),configuration.is_verbose); [recording, segment_data] = recording.get_next_window(); + % Calculating the wPLI and binarize it + segment_wpli = wpli(segment_data, number_surrogate, p_value); + result.data.wpli(i,:,:) = segment_wpli; + b_wpli = binarize_matrix(threshold_matrix(segment_wpli, threshold)); + % Calculating hub data for the segment - [channel_index, channel_degree, relative_position, graph] = hub_location(segment_data, recording.channels_location, number_surrogate, p_value, threshold); + [hub_index, weights] = binary_hub_location(b_wpli, location); % Saving the hub data for this segment - result.data.hub_index(i) = channel_index; - result.data.hub_degree(i) = channel_degree; - result.data.hub_relative_position(i) = relative_position; - result.data.graph(i,:,:) = graph; + result.data.hub_index(i) = hub_index; + result.data.hub_weights(i,:) = weights; end end diff --git a/README.md b/README.md index e6cde826..c93fd606 100644 --- a/README.md +++ b/README.md @@ -17,6 +17,21 @@ A hub is a node in the brain with a lot of incoming connection, but a few output We can use the degree to approximate this, however using betweeness centrality help a lot the analysis. Below you will see the hub location definition we've build using the [Brain Connectivity Toolbox](https://sites.google.com/site/bctnet/). +Concretely this is how we can use the `na_hub_location.m` functino: +```matlab + + recording = load_set('name_of_data.set',path_to_data); + frequency_band = [7 13]; + window_size = 10; + number_surrogate = 10; + p_value = 0.05; + threshold = 0.10; + step_size = 10; + result_hl = na_hub_location(recording, frequency_band, window_size, step_size, number_surrogate, p_value, threshold); + +``` +For more detailed help type `help na_hub_location` at the MATLAB command prompt. + The definition of a hub we are using is the following: `hub = Max of (1.0*norm_betweenes_centrality + 1.0*norm_degree)`