diff --git a/JOVE_GUI3.m b/JOVE_GUI3.m
index b30b587..3714edc 100644
--- a/JOVE_GUI3.m
+++ b/JOVE_GUI3.m
@@ -1,559 +1,559 @@
 function varargout = JOVE_GUI3(varargin)
 
 gui_Singleton = 1;
 gui_State = struct('gui_Name',       mfilename, ...
                    'gui_Singleton',  gui_Singleton, ...
                    'gui_OpeningFcn', @JOVE_GUI3_OpeningFcn, ...
                    'gui_OutputFcn',  @JOVE_GUI3_OutputFcn, ...
                    'gui_LayoutFcn',  [] , ...
                    'gui_Callback',   []);
 if nargin && ischar(varargin{1})
     gui_State.gui_Callback = str2func(varargin{1});
 end
 
 if nargout
     [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
 else
     gui_mainfcn(gui_State, varargin{:});
 end
 
 function JOVE_GUI3_OpeningFcn(hObject, eventdata, handles, varargin)
 
 handles.output = hObject;
 
 addpath(genpath(fullfile(pwd,'Utilities')));
 
 handles.isPR = 0;
 
 handles.W = realmax;
 handles.Step = realmax;
 handles.Boot = 10;
 
 handles.TC = {};
 handles.TC_sub = {};
 
 handles.CurrentState = [];
 
 handles.IsDataLoaded = 0;
 
 handles.StartTiming = 1;
 handles.EndTiming = 2;
 
 handles.ISFC = {};
 
 handles.RefSize = 2;
 
 handles.SesInfo = [];
 
 handles.SN = '';
 
 % Creates the session table (first, one row and only one value)
 set(handles.SesTable,'Data',{1});
 set(handles.SesTable,'ColumnWidth',{10});
 
 set(handles.GroupPlot,'XTick',[]);
 set(handles.GroupPlot,'YTick',[]);
 set(handles.GroupPlot,'Visible','off');
 set(handles.GroupPlot,'Box','off');
 
 guidata(hObject, handles);
 
 function varargout = JOVE_GUI3_OutputFcn(hObject, eventdata, handles) 
 varargout{1} = handles.output;
 
 
 
 function WindowText_Callback(hObject, eventdata, handles)
 
 % If the high-pass takes a reasonable value, then we validate it
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) > 1) && ...
         (str2double(get(hObject,'String')) < size(handles.TC{1},2))
     
     handles.W = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 function WindowText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter window size...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 function WindowText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 guidata(hObject,handles);
 
 
 
 
 function BootText_Callback(hObject, eventdata, handles)
 
 % If the high-pass takes a reasonable value, then we validate it
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) > 1)
     
     handles.Boot = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 % --- Executes during object creation, after setting all properties.
 function BootText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter boot folds...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 function BootText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 guidata(hObject,handles);
 
 
 
 
 
 function StepText_Callback(hObject, eventdata, handles)
 
 % If the high-pass takes a reasonable value, then we validate it
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) >= 1) && ...
         (str2double(get(hObject,'String')) < size(handles.TC{1},2))
     
     handles.Step = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 function StepText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter step size...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 function StepText_ButtonDownFcn(hObject,eventdata,handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 guidata(hObject,handles);
 
 
 
 
 function StartText_Callback(hObject, eventdata, handles)
 
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) >= 1) && ...
         (str2double(get(hObject,'String')) < size(handles.TC{1},2))
     
     handles.StartTiming = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
     
     for s = 1:length(handles.TC_sub)
         handles.TC_sub{s} = handles.TC{s}(:,handles.StartTiming:handles.EndTiming);
     end
     
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 
 guidata(hObject,handles);
 
 
 % --- Executes during object creation, after setting all properties.
 function StartText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter start timing...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 
 function StartText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 guidata(hObject,handles);
 
 
 
 function LoadButton_Callback(hObject, eventdata, handles)
 
 out_idx = [];
 
 % Selection of all the fMRI files
 [filename1,pathname1]=uigetfile({'*.*','All Files'},...
-  'Select fMRI files (NIFTI)...','MultiSelect','on');
+  'Select fMRI files...','MultiSelect','on');
 
 % If the user has indeed entered files
 if ~isequal(filename1,0) || ~isequal(pathname1,0)
     
         % The files are loaded sequentially: File{t} contains the path
         % towards the t-th frame
         for s = 1:length(filename1)   
             load((fullfile(pathname1, filename1{s})));
             assignin('base','ToSave', ToSave);
             handles.TC{s} = ToSave;
             
             [~,idx_toremove] = RemoveNaNRegions(handles.TC{s});
             
             out_idx = unique([out_idx,idx_toremove]);
         end
 
         for s = 1:length(filename1)
             
             handles.TC{s}(ismember(1:size(handles.TC{s},1),out_idx),:) = NaN;
             handles.TC_sub{s} = handles.TC{s};
         end
         
         handles.IsDataLoaded = 1;
 end
 
 handles.EndTiming = size(handles.TC{1},2);
 
 handles.CurrentState = zeros(1,length(handles.TC));
 
 n_rows = floor(length(handles.TC)/8);
 n_add = mod(length(handles.TC),8);
 
 if n_rows > 0
     
     tmp_data = num2cell(zeros(n_rows+1,8));
     
     for i = 1:n_rows
         for j = 1:8
             tmp_data{i,j} = 1;
         end
     end
 end
 
 for a = 1:n_add
     tmp_data{n_rows+1,a} = 1;
 end
 
 
 set(handles.SesTable,'Data',tmp_data);
 
 % Plots the depiction of groups
 handles.GroupPlot = JOVE_PlotISFCGroups(handles.GroupPlot,handles.Boot,handles.CurrentState);
 
 guidata(hObject,handles);
 
 
 
 
 function PlotButton_Callback(hObject, eventdata, handles)
 
 % Gets session labels
 handles.SesInfo = cell2mat(get(handles.SesTable,'Data'));
 
 tmp_ses = handles.SesInfo';
 tmp_ses = tmp_ses(:);
 tmp_ses(tmp_ses == 0) = [];
 
 if handles.isPR
     handles.TC_sub = JOVE_PhaseRandomise(handles.TC_sub);
 else
     
 end
 
 handles.ISFC = JOVE_Perform_Bootstrapping(handles.TC_sub,tmp_ses,handles.W,handles.Step,handles.TR,handles.Boot,handles.RefSize,handles.GroupPlot,handles.ISFCPlot,hObject,handles);
 
 [tmp,IDX] = sort(sum(abs(handles.ISFC{1}),2),'descend');
 IDX = IDX(~isnan(tmp));
 
 try
     handles.ISFCPlot = JOVE_PlotAtlasTC(handles.ISFCPlot,handles.ISFC{1},IDX(1:50),'ISFC',handles.W,handles.Step,handles.TR);
 catch
     errordlg('No computation for your subject of interest; run again with larger boot number !');
 end
 
 guidata(hObject,handles);
 
 
 function EndText_Callback(hObject, eventdata, handles)
 
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) >= handles.StartTiming) && ...
         (str2double(get(hObject,'String')) <= size(handles.TC{1},2))
     
     handles.EndTiming = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
     
     for s = 1:length(handles.TC_sub)
         handles.TC_sub{s} = handles.TC{s}(:,handles.StartTiming:handles.EndTiming);
     end
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 % --- Executes during object creation, after setting all properties.
 function EndText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter end timing...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 
 function EndText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 
 guidata(hObject,handles);
 
 
 
 function RefText_Callback(hObject, eventdata, handles)
 
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) > 1) && ...
         (str2double(get(hObject,'String')) < length(handles.TC))
     
     handles.RefSize = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 
 
 
 function RefText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter ref. number...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 
 
 function RefText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 guidata(hObject,handles);
 
 
 % --- Executes when entered data in editable cell(s) in SesTable.
 function SesTable_CellEditCallback(hObject, eventdata, handles)
 % hObject    handle to SesTable (see GCBO)
 % eventdata  structure with the following fields (see MATLAB.UI.CONTROL.TABLE)
 %	Indices: row and column indices of the cell(s) edited
 %	PreviousData: previous data for the cell(s) edited
 %	EditData: string(s) entered by the user
 %	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
 %	Error: error string when failed to convert EditData to appropriate value for Data
 % handles    structure with handles and user data (see GUIDATA)
 
 
 % --------------------------------------------------------------------
 function SesTable_ButtonDownFcn(hObject, eventdata, handles)
 % hObject    handle to SesTable (see GCBO)
 % eventdata  reserved - to be defined in a future version of MATLAB
 % handles    structure with handles and user data (see GUIDATA)
 
 
 % --- Executes when selected cell(s) is changed in SesTable.
 function SesTable_CellSelectionCallback(hObject, eventdata, handles)
 % hObject    handle to SesTable (see GCBO)
 % eventdata  structure with the following fields (see MATLAB.UI.CONTROL.TABLE)
 %	Indices: row and column indices of the cell(s) currently selecteds
 % handles    structure with handles and user data (see GUIDATA)
 
 
 
 function ClearButton_Callback(hObject, eventdata, handles)
 
 handles.isPR = 0;
 
 handles.W = realmax;
 handles.Step = realmax;
 handles.Boot = 10;
 
 handles.TC = {};
 handles.TC_sub = {};
 
 handles.CurrentState = [];
 
 handles.IsDataLoaded = 0;
 
 handles.StartTiming = 1;
 handles.EndTiming = 2;
 
 handles.ISFC = {};
 
 handles.RefSize = 2;
 
 handles.SesInfo = [];
 
 handles.SN = '';
 
 % Creates the session table (first, one row and only one value)
 set(handles.SesTable,'Data',{1});
 set(handles.SesTable,'ColumnWidth',{10});
 
 set(handles.GroupPlot,'XTick',[]);
 set(handles.GroupPlot,'YTick',[]);
 set(handles.GroupPlot,'Visible','off');
 set(handles.GroupPlot,'Box','off');
 
 set(handles.WindowText,'Enable','off');
 set(handles.StepText,'Enable','off');
 set(handles.TRText,'Enable','off');
 set(handles.BootText,'Enable','off');
 set(handles.RefText,'Enable','off');
 set(handles.StartText,'Enable','off');
 set(handles.EndText,'Enable','off');
 set(handles.SaveText,'Enable','off');
 
 set(handles.WindowText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.StepText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.TRText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.BootText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.RefText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.StartText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.EndText,'BackgroundColor',240/255*[1 1 1]);
 set(handles.SaveText,'BackgroundColor',240/255*[1 1 1]);
 
 WindowText_CreateFcn(handles.WindowText, eventdata, handles);
 StepText_CreateFcn(handles.StepText, eventdata, handles);
 TRText_CreateFcn(handles.TRText, eventdata, handles);
 BootText_CreateFcn(handles.BootText, eventdata, handles);
 RefText_CreateFcn(handles.RefText, eventdata, handles);
 StartText_CreateFcn(handles.StartText, eventdata, handles);
 EndText_CreateFcn(handles.EndText, eventdata, handles);
 SaveText_CreateFcn(handles.SaveText,eventdata,handles);
 
 cla(handles.GroupPlot);
 cla(handles.ISFCPlot);
 
 guidata(hObject, handles);
 
 
 
 function SaveButton_Callback(hObject, eventdata, handles)
 
 ISFC_Final = handles.ISFC;
 
 save(handles.SN,'ISFC_Final');
 
 
 
 function SaveText_Callback(hObject, eventdata, handles)
 
 if (~isempty(get(hObject,'String')))
     handles.SN = get(hObject,'String');
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 
 function SaveText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter save name...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 
 
 function SaveText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 
 guidata(hObject,handles);
 
 
 
 function TRText_Callback(hObject, eventdata, handles)
 
 % If the high-pass takes a reasonable value, then we validate it
 if (~isempty(str2double(get(hObject,'String')))) && ...
         (str2double(get(hObject,'String')) > 0)
     
     handles.TR = str2double(get(hObject,'String'));    
     set(hObject,'BackgroundColor', [0.4 0.6 0.4]);
 else 
     set(hObject,'BackgroundColor', [0.93 0.84 0.84]);
 end
 
 guidata(hObject,handles);
 
 
 function TRText_CreateFcn(hObject, eventdata, handles)
 
 set(hObject,'String','Enter TR...');
 set(hObject,'FontAngle','italic');
 
 guidata(hObject,handles);
 
 
 
 function TRText_ButtonDownFcn(hObject, eventdata, handles)
 
 set(hObject,'Enable','on');
 set(hObject,'String','');
 set(hObject,'FontAngle','normal');
 uicontrol(hObject);
 guidata(hObject,handles);
 
 
 
 function checkbox1_Callback(hObject, eventdata, handles)
 
 if get(hObject,'Value')
     handles.isPR = 1;
 else
     handles.isPR = 0;
 end
 
 uicontrol(hObject);
 guidata(hObject,handles);
diff --git a/Utilities/ISFC_computations/GetConVec.m b/Utilities/ISFC_computations/GetConVec.m
deleted file mode 100644
index 300a2ce..0000000
--- a/Utilities/ISFC_computations/GetConVec.m
+++ /dev/null
@@ -1,34 +0,0 @@
-%% This function converts a vector of regions into a vector of connections
-% for which one cell has the format 'i-j', with i and j two integer numbers
-% highlighting two atlas regions
-%
-% Inputs:
-%
-% - RV is the vector of regional indices (n_ROI x 1)
-%
-% Outputs:
-%
-% - CV is the cell array of connections (length 1 x n_conn)
-function [CV] = GetConVec(RV)
-
-    n_ROI = length(RV);
-
-    % Fills up a n_ROI x n_ROI matrix with regional indices
-    for i = 1:n_ROI
-        Mat1(i,:) = RV(i)*ones(1,n_ROI);
-    end
-
-    for i = 1:n_ROI
-        Mat2(:,i) = RV(i)*ones(n_ROI,1);
-    end
-
-    % Samples the upper diagonal matrix only
-    V1 = jUpperTriMatToVec(Mat1);
-    V2 = jUpperTriMatToVec(Mat2);
-
-    % From it, creates the output vector
-    for i = 1:length(V1)
-        CV{i} = [num2str(V1(i)),'-',num2str(V2(i))];
-    end
-
-end
\ No newline at end of file
diff --git a/Utilities/ISFC_computations/JOVE_ISFC.m b/Utilities/ISFC_computations/JOVE_ISFC.m
deleted file mode 100644
index d4ba2de..0000000
--- a/Utilities/ISFC_computations/JOVE_ISFC.m
+++ /dev/null
@@ -1,261 +0,0 @@
-%% This script runs the ISFC computations and ISFC transient extraction on
-% a set of typically developing subjects, for a movie-watching paradigm
-%
-% Code written by Thomas Bolton, last checked on July 24th 2018
-
-%% 1. Data loading, paths generation and parameters
-
-% Adds the path towards all utilities
-addpath(genpath(fullfile(pwd,'Utilities')));
-
-%%%%%%%%%%%%%%%%%% To be entered by the user %%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Name of the atlas file
-Atlas_name = 'Craddock_300.nii';
-
-% Name of the data file
-Data_name = 'JOVE_Full_Data.mat';
-
-% Name of the file containing the output data
-Save_name = 'JOVE_Results';
-
-% TR for the considered dataset
-TR=2.0;
-
-% Subjects that were included in the atlasing step a priori (i.e., name of
-% the related subject folders)
-Subjects_OI = {'S1','S9','S10','S14','S16','S17','S19','S20','S21','S23','S25',...
-    'S26','S27','S31','S32'};
-
-Ses_prefix = 'SES';
-
-% Threshold to use for motion selection (scrubbing), according to Power's
-% criterion; advised value is 0.5 [mm]
-Mot_Threshold = 0.5;
-
-% Percentage frame that can maximally be scrubbed before the subject is
-% taken out; advised value is 10 [%], can go up to 30 [%]
-TScrub = 10;
-
-% Size of the sliding windows (in TRs); advised value is W = 10 [TR]
-W = 10;
-
-% Step between two windows (in TRs); advised value is Step = 1 [TR]
-Step = 1;
-
-% Indices of the sessions, from the ones to analyze, combining RS and
-% movie-watching subparts
-combined_sessions = [1,2];
-
-% Timing (in samples from the data matrix) of the movie-watching and RS
-% subparts for those sessions
-MW_timing = 3:178;
-RS_timing = 194:340;
-
-% Number of bootstrapping folds to use for ISFC computations
-n_boot = 250;
-
-% Number of reference subjects to use for the bootstraping process
-n_refsubj = 6;
-
-% alpha-value to use for ISFC thresholding
-Alpha = 0.01;
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Path towards the atlas
-PathToAtlas = fullfile(pwd,'Utilities',Atlas_name);
-
-% Path towards the folder will all my subjects
-PathToSubjFolder = fullfile(pwd,'Subjects');
-
-% Loads the data (called 'TCS' from the previous script)
-load(fullfile(pwd,'Results',Data_name));
-
-% Number of different sessions acquired per subject
-n_sessions = size(TCS,1);
-
-% For each session, we assign a label
-Subject_Label = repmat(Subjects_OI,size(TCS,1),1);
-Subject_Label = Subject_Label(:);
-
-% Please switch to true if desiring to plot some figures of the process
-is_plot = 1;
-
-% Please switch to true if desiring to plot some text about the process
-is_text = 1;
-
-%% 2. Removal of bad quality sessions
-
-% Number of subjects
-n_subj = length(Subjects_OI);
-
-% Label for session (1 or 2)
-Session_Label = repmat(1:n_sessions,1,length(Subjects_OI));
-
-% Computation of motion for all subjects
-[~,~,PScrub] = ComputeFD(PathToSubjFolder,Subjects_OI,Mot_Threshold,Ses_prefix);
-
-% We keep subjects if they show less than 30% of time points with large
-% framewise displacement
-ToKeep = (PScrub <= TScrub);
-
-if is_text
-    disp([num2str(sum(~ToKeep)),'/',num2str(length(ToKeep)),' sessions to discard...']);
-end
-
-if is_plot
-    figure;
-    set(gcf,'color','w');
-    bar(1:length(PScrub),PScrub);
-    set(gca,'Box','off');
-    xlabel('Sessions');
-    ylabel('Percentage scrubbed frames [%]');
-    xlim([0.5 length(PScrub)+0.5]);
-    ylim([0 50]);
-    hold on;
-    plot([1 length(PScrub)],[TScrub TScrub],'k--');
-end
-
-Session_Label = Session_Label(ToKeep);
-Subject_Label = Subject_Label(ToKeep);
-
-% I convert my 3 x 30 data into a 1 x 90 format in order to keep only the
-% sessions of interest
-TCS_all = reshape(TCS,1,max(Session_Label)*size(TCS,2));
-TCS_all = TCS_all(ToKeep);
-
-clear PScrub
-
-% Removes the regions in which at least one of the subjects has NaN values
-% (i.e. the region did not exist). We use all the sessions for this step
-[TCS_ip,idx_removedROIs] = RemoveNaNRegions(TCS_all);
-
-if is_text
-   disp([num2str(length(idx_removedROIs)),' regions removed from the atlas...']); 
-end
-
-% Filtering of the time courses at cutoff 1/W [s]
-TCS_rear = cell(1,length(TCS_ip));
-for i = 1:length(TCS_ip)
-    TCS_rear{i} = y_IdealFilter(TCS_ip{i}', TR, [1/(W*TR) 0]);
-    TCS_rear{i} = TCS_rear{i}';
-end
-
-%% 3. Atlas-related details
-
-% We want a vector indexing the ROIs, going from 1 to the number of regions
-% kept (not NaN)
-ROI_vec = 1:size(TCS_all{1},1);
-ROI_vec(idx_removedROIs) = [];
-ROI_vec = 1:length(ROI_vec);
-
-% Number of regions and of connections
-n_ROI = length(ROI_vec);
-n_conn = n_ROI*(n_ROI-1)/2;
-
-if is_text
-    disp(['Total number of regions: ',num2str(n_ROI),'...']);
-    disp(['Total number of connections: ',num2str(n_conn),'...']);
-end
-
-% Converts the region vector into a connection vector (an example entry is
-% '1-13' between ROI 1 and ROI 13)
-Con_vec = GetConVec(ROI_vec);
-
-% I gather the labels of the regions involved in each connection (for the
-% above example, C1(i) = 1 and C2(i) = 13)
-[C1,C2] = ParseConVec(Con_vec);
-
-% Creation of the atlas codebook and of 3D center of gravity information,
-% with COG of size 3 x n_ROI
-[CodeBook,COG] = CreateCodeBook(PathToAtlas,ROI_vec,idx_removedROIs);
-
-% Euclidean distance between all connections considered
-for i = 1:n_conn
-    D_euclidean(i,1) = sqrt(sum((squeeze(COG(:,C1(i)))-squeeze(COG(:,C2(i)))).^2));
-end
-
-if is_plot
-    figure;
-    set(gcf,'color','w');
-    hist(D_euclidean,1000);
-    xlim([0 200]);
-    set(gca,'Box','off');
-    xlabel('Distance between regions [mm]');
-    ylabel('Number of connections [-]');
-end
-
-%% 4. Extraction of resting-state and movie-watching portions
-
-idx_RSMW = 1;
-idx_RS = 1;
-
-% We iterate through all the preprocessed data
-for i = 1:length(TCS_rear)
-    
-    % tmp contains the time courses (n_ROI x n_TPs) for one subject
-    tmp = TCS_rear{i};
-    
-    % If the current recording is a combined stimulus + resting-state case,
-    % then we sample both subparts
-    switch (ismember(Session_Label(i),combined_sessions))
-        case 0
-            TC_RS2{idx_RS} = tmp;
-            idx_RS = idx_RS + 1;
-            
-        case 1
-            TC_RS{idx_RSMW} = tmp(:,RS_timing);
-            TC_MW{idx_RSMW} = tmp(:,MW_timing);
-
-            idx_RSMW = idx_RSMW + 1;
-    end
-    
-    clear tmp
-end
-
-if is_text
-    disp(['Number of retained pure resting-state sessions: ',num2str(length(TC_RS2)),'...']);
-    disp(['Number of retained stimulus-related sessions: ',num2str(length(TC_RS)),'...']);
-end
-
-%% 5. Computing ISFC estimates on the data
-
-% First, we wish to solely perform ISFC compuations on the movie-watching
-% subparts; we create the ISFC_MW cell array, each cell about to contain
-% the outputs for one subject (n_conn x n_slidingwindow_timepoint); then,
-% we do the same for the resting-state subparts from movie sessions, and
-% finally for the 'pure resting-state' sessions
-
-if is_text
-    disp('Computing ISFC on stimulus-related sessions (stimulus segments)...');
-end
-[ISFC_MW] = Perform_Bootstrapping(TC_MW,Session_Label(ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);
-
-if is_text
-    disp('Computing ISFC on stimulus-related sessions (resting-state segments)...');
-end
-[ISFC_RS] = Perform_Bootstrapping(TC_RS,Session_Label(ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);
-
-if is_text
-    disp('Computing ISFC on purely resting-state sessions...');
-end
-[ISFC_RS2] = Perform_Bootstrapping(TC_RS2,Session_Label(~ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);
-
-%% 6. Significance thresholding of connectivity pairs
-if is_text
-    disp('Thresholding ISFC time courses...');
-end
-
-% First, I want to determine when in time, and for what connections, I
-% observe a significant excursion above from the mean connectivity
-% SignMat has size n_conn x n_TP x n_sessions, and contains:
-% -1 if there is a negative excursion
-% 0 if there is no excursion
-% +1 if there is a positive excursion
-SignMat = FindExcursions(ISFC_MW,ISFC_RS,ISFC_RS2,Alpha);
-
-% Saves the data
-save(fullfile(pwd,'Results',Save_name),'CodeBook','TC_RS','TC_MW',...
-    'TC_RS2','ISFC_MW','ISFC_RS','ISFC_RS2','SignMat','Alpha','W','Step',...
-    'n_boot','n_refsubj','TR','Mot_Threshold','TScrub','-v7.3');
\ No newline at end of file
diff --git a/Utilities/ISFC_computations/ParseConVec.m b/Utilities/ISFC_computations/ParseConVec.m
deleted file mode 100644
index e9fe246..0000000
--- a/Utilities/ISFC_computations/ParseConVec.m
+++ /dev/null
@@ -1,19 +0,0 @@
-%% This function parses a cell array of connection indices of the form
-% 'i-j' into the two integers i and j
-%
-% Inputs:
-%
-% - CV is the connection cell array (each cell has the form 'i-j')
-%
-% Outputs:
-%
-% - C1 and C2 are two n_conn long vectors containing, for each, the index
-% of one of the regions present in the (i,j) regional pair
-function [C1,C2] = ParseConVec(CV)
-
-    for i = 1:length(CV)
-        [tmp3,tmp4] = strtok(CV{i},'-');
-        C1(i) = str2num(tmp3);
-        C2(i) = str2num(tmp4(2:end));
-    end
-end
\ No newline at end of file
diff --git a/Utilities/ISFC_computations/Perform_Bootstrapping.m b/Utilities/ISFC_computations/Perform_Bootstrapping.m
deleted file mode 100644
index 3e375b2..0000000
--- a/Utilities/ISFC_computations/Perform_Bootstrapping.m
+++ /dev/null
@@ -1,129 +0,0 @@
-%% This function performs the ISFC bootstrapping scheme on a specific type
-% of session for a group of subjects
-%
-% Inputs:
-%
-% - TC contains the time courses of all sessions to jointly compute ISFC on
-% in a cell array format (each cell of size n_ROI x n_timepoints)
-% - Session_Label contains the index of the session of each recording; this
-% is used to include an equal number of samples from each session in the
-% reference group of the bootstrapping scheme
-% - W is the window length (in TR)
-% - S is the window step size (in TR)
-% - n_boot is the number of bootstrapping folds to perform
-% - n_refsubj is the number of subjects to use in the reference group at
-% each fold
-function [ISFC] = Perform_Bootstrapping(TC,Session_Label,W,S,n_boot,n_refsubj)
-
-    % Number of sliding window estimates that we will have according to W
-    % (window length) and S (step size)
-    n_swtp = floor((size(TC{1},2)-W)/S)+1;
-    
-    n_ROI = size(TC{1},1);
-    n_conn = n_ROI*(n_ROI-1)/2;
-    
-    % Session indices present in this call
-    sesvals = unique(Session_Label);
-    
-    % ISFC will contain the final ISFC estimates for all sessions 
-    % (one per cell)
-    ISFC = cell(1,length(Session_Label));
-    [ISFC{:}] = deal(zeros(n_conn,n_swtp));
-
-    % Denom_factor will highlight, for each session, how many times ISFC
-    % was computed on it (i.e., how many times it was not part of the
-    % reference group out of the n_boot folds)
-    Denom_factor = cell(1,length(Session_Label));
-    [Denom_factor{:}] = deal(0);
-
-    % tmp_denom will temporarily contain the information for a given fold,
-    % with 1 highlighting sessions on which ISFC was computed 
-    tmp_denom = cell(1,length(Session_Label));
-    [tmp_denom{:}] = deal(0);
-
-    % idx_allsessions contains the indices of all included sessions, while
-    % idx_sepsessions will contain all indices of session s in its cell s
-    idx_allsessions = 1:length(Session_Label);
-    idx_sepsessions = cell(1,length(sesvals));
-    
-    % n_recperses is the number of recordings that should be sampled, at
-    % each fold, for each session index; it can be a non-integer value at
-    % this point (e.g., 1.5 if we want to have 6 reference samples out of 4
-    % different session subtypes)
-    n_recperses = n_refsubj/length(sesvals);
-
-    % We iterate over folds...
-    for boot = 1:n_boot
-        
-        if boot == 1
-            tic;
-        end
-        
-        idx_inref = [];
-        tmp_rand = 0;
-        
-        for s = 1:length(sesvals)
-            
-            % For each session subtype s, samples only the related indices
-            % and permutes them into a random order
-            idx_sepsessions{s} = idx_allsessions(Session_Label == sesvals(s));
-            idx_sepsessions{s} = idx_sepsessions{s}(randperm(length(idx_sepsessions{s})));
-            
-            % For each session subtype s, picks how many samples should be
-            % incorporated in the reference group
-            if s == length(sesvals)
-                tmp_rand(s) = n_refsubj-sum(tmp_rand);
-            else
-                tmp_rand(s) = randi([floor(n_recperses),ceil(n_recperses)]);
-            end
-            
-            % Adds the sample indices from session subtype s to idx_inref
-            idx_inref = [idx_inref,idx_sepsessions{s}(1:tmp_rand(s))];
-        end
-        
-        clear tmp_rand
-        
-        % Sorts the selected indices in ascending order
-        idx_inref = sort(idx_inref,'ascend');
-        
-        % idx_notinref contains all the indices of the sessions on which
-        % ISFC will be computed at this fold (not in the reference)
-        idx_notinref = 1:length(TC);
-        idx_notinref(idx_inref) = [];
-
-        % Selects the sessions in the reference group, and the ones not in
-        % there
-        tmp_ref = TC(idx_inref);
-        tmp_notref = TC(idx_notinref);
-
-        % Computes ISFC for tmp_notref, using tmp_ref as a reference
-        % subgroup, for the current fold
-        ISFC_onefold = ISFC_ToReferenceGroup(tmp_notref,tmp_ref,W,S);
-        
-        % We save the current fold output in tmp_cell, and sum it to the
-        % ISFC variable (which keeps track of the sum of ISFC through all 
-        % folds)
-        tmp_cell = cell(1,length(Session_Label));
-        [tmp_cell{:}] = deal(zeros(n_conn,n_swtp));
-        tmp_cell(idx_notinref) = ISFC_onefold;
-        ISFC = cellfun(@(x,y) x+y,tmp_cell,ISFC,'UniformOutput',0);
-        clear tmp_cell
-
-        % At the same time, we also keep track of how many times each
-        % subject belonged to the non-reference side in Denom_factor, which
-        % is incremented as needed here
-        tmp_denom(idx_notinref) = {1};
-        Denom_factor = cellfun(@(x,y) x+y,tmp_denom,Denom_factor,'UniformOutput',0);
-        [tmp_denom{:}] = deal(0);
-        
-        if boot == 1
-            time_onefold = toc;
-            
-            disp(['Estimated total run time: ',num2str((n_boot-1)*time_onefold),' seconds...']);
-        end
-    end
-
-    % Averages the results; we divide by the number of times a given subject has
-    % been included
-    ISFC = cellfun(@(x,y) x/y,ISFC,Denom_factor,'UniformOutput',0);
-end
\ No newline at end of file
diff --git a/Utilities/ISFC_computations/SlidingWindow.m b/Utilities/ISFC_computations/SlidingWindow.m
deleted file mode 100755
index c7300f8..0000000
--- a/Utilities/ISFC_computations/SlidingWindow.m
+++ /dev/null
@@ -1,213 +0,0 @@
-%% This function performs sliding window analysis on the data of one 
-% subject. It can perform it the standard way, or the 'IFSC' way
-%
-% Inputs:
-% - X is the data of all the subjects that we want to consider (cell array with each cell n_ROI x n_TP)
-% - W is the window size (in samples)
-% - Step is the window Step (in samples)
-% - ConnType is the type of measure to use (e.g. 'Pearson')
-% - SWType is the type of sliding window correlation to perform
-%
-% Outputs:
-% - Y is the output (n_pairs x n_swtp)
-function [Y,CM_stat] = SlidingWindow(X,W,Step,ConnType,SWType)
-
-    % Number of time points of the data
-    n_TP = size(X{1},2);
-    
-    % Number of regions
-    n_ROI = size(X{1},1);
-    
-    % Number of subjects
-    n_subjects = length(X);
-    
-    % Number of pairs of regions
-    n_pairs = n_ROI * (n_ROI-1) / 2;
-
-    % Current end of the window; we want to iterate the computations until
-    % the window end is past the end of the data
-    WindowEnd = W;
-    WindowStart = 1;
-    
-    time_w = 1;
-    
-    full_index = 1;
-    
-    Y = cell(1,n_subjects);
-    CM_stat = cell(1,n_subjects);
-    
-    while WindowEnd <= n_TP
-        
-        %disp(['Currently in round ',num2str(full_index),' of computations']);
-        
-        switch SWType
-            % Normal FC across regions from the same subject
-            case 'FC'   
-                 
-                for i = 1:n_subjects
-                    
-                    % Samples the part of the time course of interest
-                    tmp = X{i}(:,WindowStart:WindowEnd);
-
-                    % Computes correlation
-                    switch ConnType
-                        case 'Pearson'
-                            % C is a n_ROI x n_ROI matrix
-                            C = corr(tmp');
-
-                            % Gets a vector with values
-                            Y{i}(:,time_w) = jUpperTriMatToVec(C);
-
-                        otherwise
-                            errordlg('Not yet implemented !');
-                    end    
-                    
-                    clear tmp
-                end
-  
-            % FC across the regions from different subjects
-            case 'ISFC'
-                
-                % For each subject, we want to compute the estimate for a
-                % given region pair comparing this particular region to all
-                % the other subject's paired area
-                for i = 1:n_subjects
-                    
-                    % tmp2 will contain our samples of time courses for all
-                    % subjects different from subject 1
-                    tmp2 = nan(n_ROI,W,n_subjects-1);
-                    
-                    % tmp is our data for subject i
-                    tmp = X{i}(:,WindowStart:WindowEnd);
-                    
-                    % Takes the indices of all the other subjects
-                    s_idx = 1:n_subjects;
-                    s_idx(i) = [];
-                    
-                    % Loops through those, and samples the data for them
-                    % into tmp2
-                    for j = 1:n_subjects-1
-                        tmp2(:,:,j) = X{s_idx(j)}(:,WindowStart:WindowEnd);
-                    end
-                    
-                    % Uses the right connectivity measure, and computes the
-                    % correlation between the mean time courses of all
-                    % subjects (except subject i) and subject i time
-                    % courses
-                    switch ConnType
-                            case 'Pearson'
-                                
-                                % I make my matrix symmetric
-                                tmp_corr = corr(tmp',squeeze(mean(tmp2,3))');
-                                tmp_corr_sym = (tmp_corr + tmp_corr')/2;
-                                
-                                Y{i}(:,time_w) = jUpperTriMatToVec(tmp_corr_sym); 
-                                
-                            otherwise
-                                errordlg('Not yet implemented !');
-                    end
-                    
-                    clear tmp2
-                    clear tmp
-                end    
-                
-            % Here, we want to implement the bootstrapping option from Chen
-            % et al. 2016 (NI) for a better ISFC computation
-            case 'Bootstrap'
-                
-                % Number of bootstrapping folds
-                n_boot = 100;
-
-                % We will consider each subject one after the other...
-                for i = 1:n_subjects
-                    
-                    disp(['Subject ',num2str(i),'...']);
-                    
-                    % This will contain the correlation coefficients
-                    % computed between subject i and all possible other
-                    % subjects
-                    CM_boot = zeros(n_ROI,n_ROI,n_boot);
-                
-                    % We want to run the same assessment many many times
-                    for boot = 1:n_boot
-
-                        tmp_corr_avg = zeros(n_ROI,n_ROI);
-                        
-                        % Samples the indices of the subjects to keep (with
-                        % replacement)
-                        idx_boot = randsample(n_subjects-1,n_subjects-1,true);
-
-                        % Takes the indices of all the other subjects
-                        s_idx = 1:n_subjects;
-                        s_idx(i) = [];
-
-                        % We will compute ISFC with every possible second
-                        % subject as selected from idx_boot
-                        for j = 1:length(idx_boot)
-
-                            % For the considered pair, n_ROI x n_ROI matrix of
-                            % correlation
-                            tmp_corr = corr(X{i}(:,WindowStart:WindowEnd)',X{s_idx(j)}(:,WindowStart:WindowEnd)');
-                            
-                            % I make my matrix symmetric
-                            tmp_corr_avg = tmp_corr_avg + ((tmp_corr + tmp_corr')/2); 
-                        end
-                        
-                        tmp_corr_avg = tmp_corr_avg/length(idx_boot);
-                        
-                        CM_boot(:,:,boot) = tmp_corr_avg;
-                    end
-                    
-                    % This contains the data of interest
-                    Y{i}(:,time_w) = jUpperTriMatToVec(median(CM_boot,3)); 
-                end
-        end
- 
-        % We shift the window by the step size
-        WindowStart = WindowStart + Step;
-        WindowEnd = WindowEnd + Step;
-        time_w = time_w + 1;
-        full_index = full_index + 1;
-    end
-    
-    switch SWType
-        case 'ISFC'
-            for i = 1:n_subjects
-                tmp2 = nan(n_ROI,n_TP,n_subjects-1);
-                
-                tmp = X{i};
-                
-                s_idx = 1:n_subjects;
-                s_idx(i) = [];
-                
-                for j = 1:n_subjects-1
-                    tmp2(:,:,j) = X{s_idx(j)};
-                end
-                
-                switch ConnType
-                        case 'Pearson'
-                            tmp_corr = corr(tmp',squeeze(mean(tmp2,3))');
-                            tmp_corr_sym = (tmp_corr + tmp_corr')/2;
-                            CM_stat{i} = tmp_corr_sym; 
-                        otherwise
-                            errordlg('Not yet implemented !');
-                end
-            end
-            
-            case 'FC'
-                for i = 1:n_subjects
-                    
-                    tmp = X{i};
-                    
-                    switch ConnType
-                        case 'Pearson'
-                            CM_stat{i} = corr(tmp');
-                        otherwise
-                            errordlg('Not yet implemented !');
-                    end    
-                end
-                
-        case 'Bootstrap'
-            CM_stat = [];
-    end
-end
\ No newline at end of file
diff --git a/Utilities/ISFC_thresholding/FindExcursions.m b/Utilities/ISFC_thresholding/FindExcursions.m
deleted file mode 100644
index 39b4755..0000000
--- a/Utilities/ISFC_thresholding/FindExcursions.m
+++ /dev/null
@@ -1,59 +0,0 @@
-%% This function computes a distribution of values for the considered
-% null data (for each connection pair), and determines (in
-% two-tailed fashion) the thresholds past which an excursion of
-% connectivity in the data is significant. Then, those moments are returned
-%
-% Inputs:
-%
-% - Data and Null1 are the actual and null data (cell arrays n_subject
-% with size n_pairs x n_sliding_window_timepoints) acquired during the
-% stimulus sessions
-% - Null2 contains the purely resting-state recordings, as opposed to the
-% resting-state recordings that were acquired within a stimulus-including
-% session
-% - percentile (in percent) is the level of significance
-%
-% Outputs:
-%
-% - Out is a 3D matrix of size n_conn x n_slidingwindow_timepoint x
-% subjects, filled with -1 (significant negative ISFC change), 0 (no
-% significant change), or +1 (significant positive ISFC change)
-function [Out] = FindExcursions(Data,Null1,Null2,percentile)
-
-    if(size(Data,1) ~= size(Null1,1))
-        errordlg('The null data and the actual data do not have the same dimensions !!');
-    end
-    
-    % Number of region pairs
-    n_pairs = size(Data{1},1);
-    n_subjects = length(Data);
-
-    % Aggregates the null (i.e., resting-state) data
-    Null_marginalized = [];
-    for s = 1:n_subjects
-        Null_marginalized = [Null_marginalized,Null1{s}];
-    end
-
-    n_subjects2 = length(Null2);
-
-    for s = 1:n_subjects2
-        Null_marginalized = [Null_marginalized,Null2{s}];
-    end
-
-    % Does the percentile computations
-    for s = 1:n_subjects
-        PCT{s} = prctile(Null_marginalized,[100-percentile, percentile],2);
-    end
-            
-    % Creates and fills the output, thresholded matrix Out with +-1 at the
-    % moments when there was a significant ISFC change, for a given subject
-    % and connection
-    Out = zeros(n_pairs,size(Data{s},2),n_subjects);
-    
-    for s = 1:n_subjects
-        for i = 1:n_pairs
-            Out(i,Data{s}(i,:)>PCT{s}(i,1),s) = 1;
-            Out(i,Data{s}(i,:)<PCT{s}(i,2),s) = -1;
-        end
-    end 
-end
\ No newline at end of file
diff --git a/Utilities/Plotting/Script_BrainHack.m b/Utilities/Plotting/Script_BrainHack.m
deleted file mode 100644
index 1ab518c..0000000
--- a/Utilities/Plotting/Script_BrainHack.m
+++ /dev/null
@@ -1,41 +0,0 @@
-%% Plots brain regions according to weights
-
-figure;
-set(gcf,'color','w');
-
-% Number of regions
-R = 360;
-
-% Creates the conections
-CM = zeros(R,R);
-
-% Defines the nodal sizes
-CC = w_fear;
-CC(CC < 0) = 0;
-
-T_conn = 0.1;
-
-Factor_SphereSize = 0.1;
-
-Factor_Col = 1;
-
-Exp_Sphere = 1;
-
-View = 1;
-
-Colormap_nodes = 'hot';
-Colormap_edges = 'jet';
-
-Gamma = 1;
-
-LinearWeight = 1;
-
-CA = [];
-
-PT = 'brainhack';
-
-title = 'test_plot';
-
-PlotBrainGraph_Kmeans_withtitle(gca,CM,CC,CodeBook,T_conn,...
-    Factor_SphereSize,Factor_Col,Exp_Sphere,View,Colormap_nodes,...
-    Colormap_edges,Gamma,LinearWeight,CA,PT,title)
\ No newline at end of file
diff --git a/Utilities/Plotting/Script_SF5.m b/Utilities/Plotting/Script_SF5.m
deleted file mode 100644
index d9479f7..0000000
--- a/Utilities/Plotting/Script_SF5.m
+++ /dev/null
@@ -1,214 +0,0 @@
-%% This script compares the results obtained on a restricted dataset with
-% two different null generation schemes
-
-% EXC_PR and EXC_RS
-n_subj = 12;
-
-% Distribution of the number of significant excursions across connections
-% and subjects
-
-%% Percentage of detected excursions with both null schemes
-
-for s = 1:n_subj
-    
-    % Number of counts per connection
-    Counts_RS(:,s) = sum(squeeze(SM_RS(:,:,s))~=0,2)/165*100;
-    
-    Counts_PR(:,s) = sum(squeeze(SM_PR(:,:,s))~=0,2)/165*100;
-end
-
-% Mean and STD
-Mu_RS = mean(Counts_RS,2);
-Std_RS = std(Counts_RS,[],2);
-
-Mu_PR = mean(Counts_PR,2);
-Std_PR = std(Counts_PR,[],2);
-
-
-% Representation of the counts with both null schemes
-figure;
-set(gcf,'color','w');
-plot_ci((1:44850)',[Mu_RS,Mu_RS-Std_RS,Mu_RS+Std_RS],...
-    'PatchColor', [0.4 0.6 0.8], 'PatchAlpha', 0.9, 'MainLineWidth', 1, ...
-    'MainLineColor', [0 0.4 1], 'LineWidth', 1, 'LineStyle','-', ...
-    'LineColor', 'none');
-ylim([0 30]);
-xlim([1 44850]);
-set(gca,'Box','off');
-xlabel('Connections');
-ylabel('Perc. exc. time points [%]');
-
-figure;
-set(gcf,'color','w');
-plot_ci((1:44850)',[Mu_PR,Mu_PR-Std_PR,Mu_PR+Std_PR],...
-    'PatchColor', [0.8 0.4 0], 'PatchAlpha', 0.9, 'MainLineWidth', 1, ...
-    'MainLineColor', [1 0.2 0.2], 'LineWidth', 1, 'LineStyle','-', ...
-    'LineColor', 'none');
-ylim([0 30]);
-xlim([1 44850]);
-set(gca,'Box','off');
-xlabel('Connections');
-ylabel('Perc. exc. time points [%]');
-
-
-% Representation of the counts with both null schemes
-figure;
-set(gcf,'color','w');
-plot_ci((1:44850)',[Mu_RS,Mu_RS-Std_RS,Mu_RS+Std_RS],...
-    'PatchColor', [0.4 0.6 0.8], 'PatchAlpha', 0.9, 'MainLineWidth', 1, ...
-    'MainLineColor', [0 0.4 1], 'LineWidth', 1, 'LineStyle','-', ...
-    'LineColor', 'none');
-ylim([0 30]);
-xlim([1 897]);
-set(gca,'Box','off');
-xlabel('Connections');
-ylabel('Perc. exc. time points [%]');
-
-
-figure;
-set(gcf,'color','w');
-plot_ci((1:44850)',[Mu_PR,Mu_PR-Std_PR,Mu_PR+Std_PR],...
-    'PatchColor', [0.8 0.4 0], 'PatchAlpha', 0.9, 'MainLineWidth', 1, ...
-    'MainLineColor', [1 0.2 0.2], 'LineWidth', 1, 'LineStyle','-', ...
-    'LineColor', 'none');
-ylim([0 30]);
-xlim([1 897]);
-set(gca,'Box','off');
-xlabel('Connections');
-ylabel('Perc. exc. time points [%]');
-
-
-%% "Distance" between the outcomes from both schemes
-
-%%%%%% Distance counting all data points
-
-% % Distance between schemes
-% Dis = abs(SM_RS) - abs(SM_PR);
-% Dis2 = squeeze(sum(Dis,2)/165);
-% 
-% % Representation of the counts with both null schemes
-% figure;
-% plot_ci((1:44850)',[mean(Dis2,2),mean(Dis2,2)-std(Dis2,[],2),mean(Dis2,2)+std(Dis2,[],2)],...
-%     'PatchColor', [0.8 0.6 0.8], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ...
-%     'MainLineColor', [0.6 0 0.6], 'LineWidth', 1, 'LineStyle','-', ...
-%     'LineColor', 'none');
-% ylim([0 0.5]);
-% 
-% mean(mean(Dis2))
-% std(mean(Dis2,2))
-
-
-%%%%%% Distance counting only data points where at least one scheme gives a
-%%%%%% significant change
-
-
-% tmp_PR = squeeze(sum(abs(SM_PR),3))/n_subj;
-% tmp_RS = squeeze(sum(abs(SM_RS),3))/n_subj;
-% 
-% Idx_PR = (tmp_PR > 0);
-% Idx_RS = (tmp_RS > 0);
-% 
-% Idx = Idx_PR | Idx_RS;
-% 
-% 
-% for s = 1:n_subj
-%    tmp1 = squeeze(SM_RS(:,:,s));
-%    tmp1 = tmp1(:);
-%    tmp1 = tmp1(Idx(:));
-% 
-%    tmp2 = squeeze(SM_PR(:,:,s));
-%    tmp2 = tmp2(:);
-%    tmp2 = tmp2(Idx(:));
-%     
-%    Dis_spec(s) = sum(abs(tmp1-tmp2))/length(tmp1);
-% end
-
-
-
-
-
-
-%% Distance for each subject, for time points across connections detected 
-% by at least one of the schemes
-for s = 1:n_subj
-
-   tmp1 = squeeze(SM_RS(:,:,s));
-   
-   Idx_RS = (abs(tmp1) > 0); 
-   
-
-   tmp2 = squeeze(SM_PR(:,:,s));
-   
-   Idx_PR = (abs(tmp2) > 0); 
-   
-   
-   Idx = Idx_RS | Idx_PR;
-   
-   tmp1 = tmp1(:);
-   tmp1 = tmp1(Idx(:));
-   
-   tmp2 = tmp2(:);
-   tmp2 = tmp2(Idx(:));
-    
-   Dis_spec(s) = sum(abs(tmp1)-abs(tmp2))/length(tmp1);
-end
-
-
-
-%% Distance between subjects (RS scheme)
-
-% n_subjpairs = n_subj* (n_subj-1) / 2;
-% 
-% pair = 1;
-% 
-% for s1 = 1:n_subj
-%     for s2 = 1:n_subj
-%         if s2 > s1
-%             
-%             tmp1 = abs(squeeze(SM_RS(:,:,s1)));
-%             tmp1b = abs(squeeze(SM_RS(:,:,s2)));
-%             
-%             tmp2 = abs(squeeze(SM_PR(:,:,s1)));
-%             tmp2b = abs(squeeze(SM_PR(:,:,s2)));
-%             
-%             IdxRS = (tmp1 > 0);
-%             IdxRSb = (tmp1b > 0);
-%             
-%             IdxPR = (tmp2 > 0);
-%             IdxPRb = (tmp2b > 0);
-%             
-%             IdxRS_final = IdxRS | IdxRSb;
-%             IdxRS_final = IdxRS_final(:);
-%             
-%             IdxPR_final = IdxPR | IdxPRb;
-%             IdxPR_final = IdxPR_final(:);
-%             
-%             tmp1 = tmp1(:);
-%             tmp2 = tmp2(:);
-%             
-%             tmp1b = tmp1b(:);
-%             tmp2b = tmp2b(:);
-%            
-%             Dis_RS(pair) = sum(abs(tmp1(IdxRS_final)-tmp1b(IdxRS_final)))/sum(IdxRS_final);
-%             Dis_PR(pair) = sum(abs(tmp2(IdxPR_final)-tmp2b(IdxPR_final)))/sum(IdxPR_final);
-% 
-%             pair = pair + 1;
-%         end
-%     end
-% end
-% 
-% % Representation of the counts with both null schemes
-% figure;
-% plot_ci((1:44850)',[mean(Dis_RS,2),mean(Dis_RS,2)-std(Dis_RS,[],2),mean(Dis_RS,2)+std(Dis_RS,[],2)],...
-%     'PatchColor', [0.4 0.6 0.8], 'PatchAlpha', 0.9, 'MainLineWidth', 1, ...
-%     'MainLineColor', [0 0.4 1], 'LineWidth', 1, 'LineStyle','-', ...
-%     'LineColor', 'none');
-% ylim([0 0.5]);
-% 
-% % Representation of the counts with both null schemes
-% figure;
-% plot_ci((1:44850)',[mean(Dis_PR,2),mean(Dis_PR,2)-std(Dis_PR,[],2),mean(Dis_PR,2)+std(Dis_PR,[],2)],...
-%     'PatchColor', [0.8 0.4 0], 'PatchAlpha', 0.9, 'MainLineWidth', 1, ...
-%     'MainLineColor', [1 0.2 0.2], 'LineWidth', 1, 'LineStyle','-', ...
-%     'LineColor', 'none');
-% ylim([0 0.5]);
\ No newline at end of file
diff --git a/Utilities/Preprocessing/ComputeFD.m b/Utilities/Preprocessing/ComputeFD.m
deleted file mode 100755
index 949b830..0000000
--- a/Utilities/Preprocessing/ComputeFD.m
+++ /dev/null
@@ -1,65 +0,0 @@
-%% This function computes framewise displacement from an rp file
-function [FD,BinMat,PScrub,MOT] = ComputeFD(PathToSubjFolder,SOI,MotThresh,Ses_prefix)
-
-    cd(PathToSubjFolder);
-    
-    idx = 1;
-    
-    for s = 1:length(SOI)
-        
-        cd(SOI{s});
-        
-        tmp_ses = dir([Ses_prefix,'*']);
-
-        for ses = 1:3
-            
-            try
-                cd(tmp_ses(ses).name);
-                
-                % Gets the motion information (two motion file names). We want to
-                % get the data from the first one
-                tmp_mot = cellstr(spm_select('List',pwd,['^' 'rp_f' '.*\.' 'txt' '$']));
-
-                % Mot is a 340 x 6 matrix containing the info the the session
-                % considered for all 6 motion parameters
-                Mot = textread(tmp_mot{1});
-                Mot = Mot(:,1:6);
-                
-               % Converts last columns into mm
-                Mot(:,4:6) = Mot(:,4:6) * 50;
-
-                % Gets the differential
-                DMot = diff(Mot,1);
-                DMot = [zeros(1,6);DMot];
-
-                % framewise displacement for the considered sessions
-                FD{idx} = sum(abs(DMot),2);
-                MOT{idx} = Mot;
-                    
-                cd('..');
-                idx = idx + 1;
-            catch
-                disp('Could not find third session...');
-                
-                FD{idx} = NaN;
-                MOT{idx} = NaN;
-                
-                idx = idx + 1;
-            end
-        end 
-        
-        cd('..');
-    end
-    
-    for i = 1:length(FD)
-        if ~isnan(FD{i})
-            PScrub(i) = 100*(sum((FD{i} > MotThresh)))/(length(FD{i}));
-        else
-            PScrub(i) = NaN;
-        end
-    end
-    % Summarizing matrix for frames to scrub 
-    BinMat = 0;
-    
-    cd('..');
-end
\ No newline at end of file
diff --git a/Utilities/Preprocessing/JOVE_Preprocessing.m b/Utilities/Preprocessing/JOVE_Preprocessing.m
deleted file mode 100644
index 6623b0f..0000000
--- a/Utilities/Preprocessing/JOVE_Preprocessing.m
+++ /dev/null
@@ -1,353 +0,0 @@
-%% This script contains the main preprocessing steps required to run the 
-% ISFC transient excursion detection methodology. In more details, the
-% following steps are applied to the data:
-%
-% - Detrending of the data
-% - Regressing out covariates (constant,linear, quadratic terms, average
-% CSF and average WM signals)
-% - Atlasing into selected atlas
-% - Scrubbing
-%
-% 
-%
-%
-%
-% Written by Thomas Bolton, last modified on July 30th 2018
-
-%% Parameters to be changed by the user
-
-%%%%%%%%% Move to before the "Subjects" folder
-
-% Name of all the subject folders through which to operate
-Subjects_OI = {'S17'};
-
-% Prefix of the folders containing the different sessions for a subject
-SES_Prefix = 'SES';
-
-% Name of the folder containing all the subjects folders
-SFold_Name = 'Subjects';
-
-% End of the file name of one functional file for a subject
-Func_Suffix = '001-01.hdr';
-
-Func_Prefix = 'rf';
-
-% File type (functional files), between 'nii' and 'hdr'
-Func_Type = 'hdr';
-
-% Prefix of the folder containing the structural data
-Struct_Prefix = 'S';
-
-Atlas_Prefix = 'wCraddock_300';
-
-Save_Name = 'JOVE_Example';
-
-% Framewise displacement threshold past which data is considered corrupted
-% (according to Power's criterion)
-Mot_Threshold = 0.5;
-
-%% Start of the preprocessing code
-
-% Turn to true if you wish to observe different outputs from the script
-is_plot = 1;
-
-% Turn to true if you wish to display different text summaries of the
-% preprocessing steps
-is_text = 1;
-
-% Current folder (at the start of the folder tree) is added with all
-% subfolders
-Now = pwd;
-addpath(genpath(pwd));
-
-% Gets to the subjects-containing folder
-cd(SFold_Name);
-
-% Iterates across subjects...
-for i = 1:length(Subjects_OI)
-    
-    if is_text
-        disp(['Preprocessing subject ',num2str(i),'/',num2str(length(Subjects_OI)),'...']);
-    end
-    
-    % Moves to the subject folder
-    cd(Subjects_OI{i});
-    
-    % Looks for the directories of interest containing the different
-    % sessions to look at (all start with 'S')
-    DOR = dir([SES_Prefix,'*']);
-    
-    % Iterates across sessions
-    for j = 1:length(DOR)
-        
-        if is_text
-            disp(['Processing session ',num2str(j),'/',num2str(length(DOR)),'...']);
-        end
-        
-        % Gets to the session directory of interest
-        cd(DOR(j).name);
-        
-        % Locates one functional file from the session in question
-        tmp_refname = dir(['*',Func_Suffix]);
-        
-        % Gets the full path of that file
-        REF = fullfile(pwd,tmp_refname(1).name);
-        
-        % Looks for the directory containing structural information and
-        % segmentation results
-        SDIR = dir([Struct_Prefix,'*']);
-        
-        % Goes to the structural directory
-        cd(SDIR(1).name);
-        
-        % Gets the link towards the native space Craddock atlas
-        Atlas_link = cellstr(spm_select('List',pwd,['^', Atlas_Prefix, '.*\.', 'nii', '$']));
-        
-        % Makes this link a full path
-        Atlas_link = fullfile(pwd,Atlas_link);
-        
-        if is_text
-           disp('Creating WM and CSF masks...'); 
-        end
-        
-        % Gets the links towards the native space WM and CSF probabilistic
-        % maps (given the c2 and c3 prefixes by SPM12)
-        WM_link = cellstr(spm_select('List',pwd,['^' 'c2' '.*\.' 'nii' '$']));
-        CSF_link = cellstr(spm_select('List',pwd,['^' 'c3' '.*\.' 'nii' '$']));
-        
-        % Converts the resolution to be of right size
-        WM = mapVolumeToVolume(WM_link{1},REF);
-        CSF = mapVolumeToVolume(CSF_link{1},REF);
-        
-        % Binarises the two volumes in very conservative manner
-        WM(WM < 0.99) = 0;
-        WM(WM >= 0.99) = 1;
-        WM = logical(WM);
-        
-        CSF(CSF < 0.99) = 0;
-        CSF(CSF >= 0.99) = 1;
-        CSF = logical(CSF);
-        
-        if is_text
-            disp(['There are ',num2str(sum(sum(sum(WM)))), ' WM voxels',...
-                ' and ',num2str(sum(sum(sum(CSF)))),' CSF voxels...']);
-        end
-        
-        % Now goes to the functional data folder to perform the operations
-        % themselves, i.e. detrending and regression (and atlasing)
-        
-        cd('..');
-
-        % Gets the links of all functional files
-        FLinks = dir([Func_Prefix,'*',Func_Type]);
-        
-        % Number of frames
-        n_frames = length(FLinks);
-        
-        if is_text
-            disp(['There are ',num2str(n_frames),' frames...']);
-        end
-        
-        if is_text
-           disp('Constructing data matrix...'); 
-        end
-        
-        % Will now read the data for that session frame by frame, and
-        % create the WM and CSF regressors
-        for d = 1:n_frames
-    
-            tmp_V = spm_read_vols(spm_vol(FLinks(d).name));
-            
-            % Regressors for WM and CSF are constructed (average of the
-            % retained voxels)
-            CSF_cov(d) = mean(tmp_V(CSF));
-            WM_cov(d) = mean(tmp_V(WM));
-            
-            % Data matrix construction (n_tp x n_vox)
-            Data(d,:) = tmp_V(:);
-            
-            clear tmp_V
-        end
-        
-        % Demeaning of regressors
-        CSF_cov = CSF_cov - mean(CSF_cov);
-        WM_cov = WM_cov - mean(WM_cov);
-
-        %% Detrending
-        
-        if is_text
-           disp('Detrending...'); 
-        end
-        
-        % Data detrending
-        Data = detrend(Data);
-        
-        % Construction of the matrix of regressors: constant, linear,
-        % quadratic, CSF and WM (no motion because scrubbing is done
-        % afterwards)
-        X = [ones(n_frames,1),(1:n_frames)',((1:n_frames).^2)',CSF_cov',WM_cov'];
-        
-        % Will contain the preprocessed data
-        Data_preproc = nan(size(Data));
-        
-        %% Regressing out covariates
-        
-        if is_text
-           disp('Regressing out covariates...'); 
-        end
-        
-        % Regresses voxel after voxel (DPARSFA toolbox function) with also
-        % adding back the mean as in the lab scripts
-        for vox = 1:size(Data,2)
-            [~,Data_preproc(:,vox)] = y_regress_ss(Data(:,vox),X);
-            Data_preproc(:,vox) = Data_preproc(:,vox) + mean(Data(:,vox));
-        end
-            
-        clear Data
-        clear X
-        
-        %% Atlasing
-        
-        if is_text
-            disp('Atlasing...');
-        end
-                
-        % Converts atlas into right resolution, and makes it 1D
-        ATLAS = mapVolumeToVolume(Atlas_link{1},REF);
-        
-        ATLAS = ATLAS(:);
-        ATLAS(isnan(ATLAS)) = 0;
-
-        % Atlases the data
-        [TC,n_disc_TC,n_kept_TC] = Make_Atlasing(Data_preproc',ATLAS);
-
-        clear ATLAS
-    
-        %% Scrubbing
-        
-        if is_text
-            disp('Scrubbing...');
-        end
-        
-        % I now want to read motion parameters
-        motfile_name = dir(['rp','*','.txt']);
-        motfile_name = motfile_name(1).name;
-        
-        % My 6 motion parameters
-        Mot = textread(motfile_name);
-        Mot = Mot(:,1:6);
-        
-        % Converts the rotational components into [mm]
-        Mot(:,4:6) = 50*Mot(:,4:6);
-        
-        % Computes FD
-        FD = sum(abs([0 0 0 0 0 0; diff(Mot)]),2);
-        
-        % Locates the time points to scrub (set to 1)
-        Scrub_mask1 = logical(FD > Mot_Threshold);
-        
-        % Locates the time point that immediately follows as well
-        Scrub_mask2 = circshift(Scrub_mask1,[1,0]);
-        Scrub_mask2(1) = 0;
-        
-        % Makes an OR operations (sets Scrub_mask element to 1 as soon as
-        % one of Scrub_mask1/Scrub_mask2 is equal to 1)
-        Scrub_mask = Scrub_mask1 | Scrub_mask2;
-        
-        if is_text
-            disp(['There are ',num2str(100*sum(Scrub_mask)/n_frames),...
-                ' frames to srub...']);
-        end
-        
-        % Finds the values at the time points that we know are not 
-        % corrupted
-        TCon = TC(:,find(~Scrub_mask));
-        TCon(isnan(TCon)) = 0;
-        
-        % tinter is all the time points
-        tinter = 1:length(Scrub_mask);
-
-        % torig is the time points that we know
-        torig = tinter(find(~Scrub_mask));
-
-        % We interpolate all time point values using the info that we
-        % know (TCon at torig)
-        TC2 = interp1(torig,TCon',tinter,'spline')';
-        
-        % We want to keep the 'NaN's that we had because they are used to
-        % remove bad regional time courses later on
-        TC2(isnan(TC)) = NaN;
-        
-        % Plotting of different aspects of the preprocessing
-        if is_plot
-            
-            % Plot of WM/CSF time courses that were regressed out, and
-            % motion parameters
-            figure;
-            set(gcf,'color','w');
-            subplot(2,1,1);
-            plot(Mot(:,1),'color',[0 0 0.6],'LineWidth',2);
-            hold on;
-            title(['Subject ',num2str(i),', session ',num2str(j)]);
-            plot(Mot(:,2),'color',[0 0.4 1],'LineWidth',2);
-            plot(Mot(:,3),'color',[0.6 0.6 1],'LineWidth',2);
-            plot(Mot(:,4),'color',[0.8 0 0.2],'LineWidth',2);
-            plot(Mot(:,5),'color',[1 0.2 0.2],'LineWidth',2);
-            plot(Mot(:,6),'color',[1 0.6 0.2],'LineWidth',2);
-            legend('X','Y','Z','Alpha','Beta','Gamma');
-            xlabel('Time');
-            ylabel('Motion [mm]');
-            set(gca,'Box','off');
-            
-            subplot(2,1,2);
-            plot(CSF_cov,'color',[0.2 0.6 0],'LineWidth',2);
-            hold on;
-            plot(WM_cov,'color',[0.6 0.6 0.6],'LineWidth',2);
-            legend('CSF','WM');
-            xlabel('Time');
-            ylabel('Intensity [a.u.]');
-            set(gca,'Box','off');    
-            
-            figure;
-            set(gcf,'color','w');
-            subplot(2,1,1);
-            hold on;
-            set(gca,'Box','off');
-            bar(n_disc_TC);
-            xlabel('Regions');
-            ylabel('Number of discarded time courses');
-            subplot(2,1,2);
-            hold on;
-            set(gca,'Box','off');
-            bar(n_kept_TC);
-            xlabel('Regions');
-            ylabel('Number of kept time courses');
-        end
-        
-        clear Mot
-        clear torig
-        clear tinter
-        clear TCon
-        clear Scrub_mask
-        clear FD
-        clear CSF_cov
-        clear WM_cov     
-        clear n_frames
-        
-        % Saves the outputs as cell arrays of time courses
-        TCS{j,i} = TC2;
-
-        cd('..');
-        
-    end
-    
-    cd('..');
-    
-end
-
-cd('..');
-
-% Saves the data in the 'Results' folder
-save(fullfile(pwd,'Results',Save_Name),'TCS');
-
-
diff --git a/Utilities/Preprocessing/RemoveNaNRegions.m b/Utilities/Preprocessing/RemoveNaNRegions.m
deleted file mode 100644
index 469095a..0000000
--- a/Utilities/Preprocessing/RemoveNaNRegions.m
+++ /dev/null
@@ -1,43 +0,0 @@
-%% This function verifies whether a given atlas time course is acceptable
-% or not
-%
-% Inputs:
-%
-% - TCin is the input cell array containing the data of each subject in a
-% cell (of size n_ROI x n_timepoints)
-%
-% Outputs:
-%
-% - TCout is the same array, having removed the bad regional time courses
-% - idx_toremove summarizes which regional indices have been removed
-%
-% Written by Thomas Bolton, last checked on July 24th
-function [TCout,idx_toremove] = RemoveNaNRegions(TCin)
-
-    % Indices of the ROIs where in at least one subject, there are NaN
-    % values
-    idx_toremove = [];
-    
-    n_ROI = size(TCin,1);
-    
-    % If there is at least a NaN value in a given time course, we include
-    % it in the ones to remove
-    if ~isempty(TCin)
-        for j = 1:n_ROI
-            if(sum(isnan(TCin(j,:))) > 0)
-                idx_toremove = [idx_toremove,j];
-            end
-        end
-    end
-    
-    % Makes each ROI appear only once in the vector
-    idx_toremove = unique(idx_toremove);
-    
-    % Only keeps the time courses that were deemed OK by the check
-    if ~isempty(TCin)
-        TCout = TCin(~ismember(1:n_ROI,idx_toremove),:);
-    else
-        TCout = NaN;
-    end
-    
-end
\ No newline at end of file
diff --git a/Utilities/Preprocessing/mapVolumeToVolume2.m b/Utilities/Preprocessing/mapVolumeToVolume2.m
deleted file mode 100644
index 0efc3b6..0000000
--- a/Utilities/Preprocessing/mapVolumeToVolume2.m
+++ /dev/null
@@ -1,45 +0,0 @@
-function Vout=mapVolumeToVolume2(Vin_fn,Vout_i)
-% map voxels in the space of input volume to voxels in the space
-% of output volumes
-% Vout_fn: uses only header info
-%
-% v1.0 Jonas Richiardi
-% - initial release, based on code by Dimitri Van De Ville and Jonas
-% Richiardi
-
-% Read output space header
-% Vout_i=spm_vol(Vout_fn);
-Vout=zeros(Vout_i.dim);
-
-% read and load input space header and data
-Vin_i=spm_vol(Vin_fn);
-Vin=spm_read_vols(Vin_i);
-
-% generate all coordinates in output space
-[x1,x2,x3]=ndgrid(1:Vout_i.dim(1),1:Vout_i.dim(2),1:Vout_i.dim(3));
-idx=1:numel(Vout); % map all voxels
-
-% take every voxel in the volume spanned by the output images,
-% compute its real-world position in mm, then map input image
-
-oobList=zeros(0,4); % list of out-of-bound input voxels
-for iter=1:length(idx),
-    oob=false;
-    % recover world-space position of this voxel in mm from affine
-    % transform matrix
-    mm=Vout_i.mat*[x1(idx(iter)) x2(idx(iter)) x3(idx(iter)) 1]';
-    % convert this position into index of the closest structural voxel
-    vx=round(Vin_i.mat\[mm(1) mm(2) mm(3) 1]'); 
-    vx(vx<=0)=1;
-    vxOri=vx;
-    % remap out-of-bounds voxels to last  
-    if vx(1)>Vin_i.dim(1), vx(1)=Vin_i.dim(1); oob=true; end 
-    if vx(2)>Vin_i.dim(2), vx(2)=Vin_i.dim(2); oob=true; end
-    if vx(3)>Vin_i.dim(3), vx(3)=Vin_i.dim(3); oob=true; end
-    if (oob==true), oobList(end+1,:)=vxOri; end
-    % idx(iter): current voxel
-    Vout(idx(iter))=Vin(vx(1),vx(2),vx(3));
-    if any(Vout(idx(iter))<0)  %Vout
-        warning('mapV2V:negativeVal',['Negative voxel values at ' num2str(iter)]);
-    end
-end;
\ No newline at end of file