%% This script runs a sliding window analysis on resting-state fMRI data,
% and subsequently applies eigenconnectivity analysis to retrieve
% connectivity building blocks
%
% Authors: Anjali Tarun, Thomas Bolton
%
% Created for the BrainHack symposium in Zurich, Friday March 3rd 2017
addpath(genpath('Functions'));
addpath(genpath('StartData'));
addpath(genpath('Utilities'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 1. Data loading and parameters %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%% A. Loading
% SW_data is a cell array, with each cell one subject data (n_ROI x
% n_timepoints)
load SW_data;
% CB is a n_ROI cell array with the short names of atlas areas; used for
% plotting
load CB;
%%% B. Parameters to leave unchanged
% Number of brain regions contained in the atlas (AAL without Pallidum)
n_ROI = 88;
% Number of subjects
n_subjects = length(SW_data);
% Number of time points
n_timepoints = size(SW_data{1},2);
% Number of eigenconnectivity building blocks to compute
n_components = 20;
%%% C. Parameters that can be changed
% Window length (in TRs)
window_length = 60;
% Window step (in TRs)
window_step = 1;
% Type of window shape desired (0: rectangular, 1: tukey, 2: gaussian,
% 3: exponential)
window = 0;
% Density of retained edges for dynamic graph analysis
rho = 10/100;
% Perform or not the demeaning for eigenconnectivity generation
% (1 = with, 0 = without)
is_demean = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% 2. Computation of connectivity time courses %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% First time span (for first window)
t_start = 1;
t_end = window_length;
% Will contain connectivity time courses for all subjects
CM = cell(1,n_subjects);
% Running across subjects...
for s = 1:n_subjects
disp(['Currently considering subject ',num2str(s),...
' for connectivity computations...']);
% We compute connectivity as long as we can still move the temporal
% subwindow forward...
ind_window = 1;
while t_end <= n_timepoints
% Data from the presently considered temporal window
tmp_data = SW_data{s}(:,t_start:t_end);
% Computation of connectivity for all region pairs according to the
% chosen method
switch window
% Rectangular window (uniform weights of 1)
case 0
Window_square = ones(1,window_length);
CM{s}(:,:,ind_window) = weightedcorrs(tmp_data',...
Window_square');
% Tukey window
case 1
r=0.5;
Window_tukey = tukeywin(window_length,r);
Window_tukey(Window_tukey==0) = realmin;
CM{s}(:,:,ind_window) = weightedcorrs(tmp_data', ...
Window_tukey);
% Gaussian window
case 2
t = linspace(0,window_length,window_length);
sigma = window_length/2;
Window_gaussian = normpdf(t,window_length/2,sigma);
Window_gaussian(Window_gaussian==0) = realmin;
CM{s}(:,:,ind_window) = weightedcorrs(tmp_data',...
Window_gaussian');
% Exponential window
case 3
t = linspace(0,window_length, window_length);
theta = window_length/3;
Window_exp = ((1-exp(-1/theta))/...
(1-exp(-window_length/theta)))*...
exp((t-window_length)/theta);
Window_exp(Window_exp==0) = realmin;
CM{s}(:,:,ind_window) = weightedcorrs(tmp_data',...
Window_exp');
end
% Shift of the temporal window by a step size
t_end = t_end + window_step;
t_start = t_start + window_step;
ind_window = ind_window + 1;
end
% After each subject is processed, we reset indices for the next one
t_end = window_length;
t_start = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% 3. Dynamic graph analysis and eigenconnectivities %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Will contain the concatenated data to run eigenconnectivity analysis on
Data = [];
% For each subject...
for s = 1:n_subjects
disp(['Currently considering subject ',num2str(s),...
' for graph/eigenconnectivity computations...']);
% Will contain the concatenated data from each subject
Data_s = [];
% For each sliding window...
for sw = 1:size(CM{s},3)
% Temporarily evaluated window (i.e. information from one connectivity matrix at one time point)
SW = squeeze(CM{s}(:,:,sw));
%%% A. Computation of graph metrics
% Clustering coefficient (computed on binarized matrix with edge density rho)
CC(s,sw,:) = ComputeC(BinarizeAdjacency(SW,rho));
% Degree (computed on binarized matrix with edge density rho)
k(s,sw,:) = sum(BinarizeAdjacency(SW,rho),1);
% Construction of the subject data matrix; Data_s has size n_conpairs x n_slidingwindowspersubject
Data_s = [Data_s,jUpperTriMatToVec(SW)];
end
% If we want to perform the subject-wise demeaning, it is applied here
if is_demean
Data_s = Data_s - repmat(mean(Data_s,2),1,size(CM{s},3));
end
% Construction of the full data matrix, of size n_conpairs x n_slidingwindowtimessubjects
Data = [Data,Data_s];
end
%%% B. Computation of eigenconnectivities
% Eigenconnectivities are simply the principal directions of the data
% matrix, contained in u
[u,S,v]=svd(Data,'econ');
% SS is a vector containing the singular values
SS = diag(S);
% If needed, we switch a given eigenconnectivity pattern so that its
% maximal value is positive. Because they index a principal direction,
% the sign of the patterns is arbitrary (it combines with the sign of the
% weights, stored in v)
for iter=1:size(u,2),
[~,idx]=max(abs(u(:,iter)));
if sign(u(idx,iter))<0,
u(:,iter)=-u(:,iter);
v(:,iter)=-v(:,iter);
end;
end;
% Plots all eigenconnectivities
myVisEigenconn(u(:,1:n_components),n_ROI);
% Example plotting command to plot one eigenconnectivity
%View_A(jVecToUpperTriMat(u(:,1),n_ROI),n_ROI,CB);