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

@@ -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

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)

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
 

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

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
 

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)`