% For reproducibility, create a file named "mat" in the workspace
% containing the same values as the Excel file provided in the
% Supplementary Materials

factor_names = {
    'Δ TOL-R', 'EDA stress induction','percent baseline',...
    'percent stress induction','HR at baseline','HR stress induction',...
    'HR hypnosis','HR stress induction 2','Scr baseline'
};

% Symmetric correlation matrix (undirected network)
cor_matrix = mat + mat';
cor_matrix(logical(eye(size(cor_matrix)))) = 1;

% Threshold
threshold = 0.4;
adj_matrix = abs(cor_matrix) > threshold;

% Create the undirected graph
G = graph(adj_matrix, factor_names,'upper');

% Number of connected components
num_components = max(conncomp(G));
disp(['Number of connected components: ', num2str(num_components)]);
% Calculate betweenness centrality
betweenness_centrality = centrality(G, 'betweenness');

disp('Betweenness centrality of each node:');
disp(betweenness_centrality);

% Indices of central nodes (betweenness > 0)
idx_cent = find(betweenness_centrality > 0);
figure;
p = plot(G, 'NodeLabel', G.Nodes.Name);
title('Network physiology of performance');

% Highlight central nodes in red
highlight(p, idx_cent,'NodeColor','r','MarkerSize',9,'LineWidth',1.5);

% GLOBAL NETWORK METRICS

% Geodesic distance matrix
D = distances(G);

% Remove Inf values (disconnected nodes) for averaging
D_noInf = D(~isinf(D));

% Compute the mean shortest path length between all node pairs
avg_distance = mean(D_noInf);
disp(['Average distance between nodes: ', num2str(avg_distance)]);

% Compute the network diameter as the maximum finite shortest path length

diameter_value = max(D_noInf);
disp(['Network diameter: ', num2str(diameter_value)]);

% Compute the network density
n = numnodes(G);
E = numedges(G);
density_value = E / (n * (n - 1));
disp(['Density: ', num2str(density_value)]);

% Closeness centrality

closeness_centrality = centrality(G, 'closeness');
disp('Closeness centrality of each node:');
disp(closeness_centrality);

