diff --git a/Code/GSP_Script_ISBI_2017_Simulated.m b/Code/GSP_Script_ISBI_2017_Simulated.m
new file mode 100644
index 0000000..55b0079
--- /dev/null
+++ b/Code/GSP_Script_ISBI_2017_Simulated.m
@@ -0,0 +1,238 @@
+%% 1. Loading of required data
+
+addpath(genpath('Utilities'));
+
+cd('..');
+cd('Data/100307');
+load('X');
+load('W');
+cd('..');
+cd('..');
+cd('Code');
+
+% Defines the HRF
+A.dt = 0.72;
+A.name = 'hrf';
+HRF = spm_get_bf(A);
+
+PathToReg = '/Volumes/Lorena/GSP_HCP_Test_Subject/100307/Task/Events';
+
+[time,Regressor] = GSP_Paradigm(PathToReg,10);
+Regressor(Regressor == 5) = 1;
+Regressor(Regressor == 6) = 2;
+Regressor(Regressor == 7) = 3;
+Regressor(Regressor == 8) = 4;
+
+% Convolves each non-constant regressor with the HRF
+U.u = Regressor';
+U.name = {'Whatever'}; % U needs a name, but the string has no meaning
+
+Regressor_real = spm_Volterra(U, HRF.bf);
+
+
+%% 2. Analysis at the level of independent Slepian vectors
+
+% To find the interesting nodes, I do a GLM beforehand
+[SIGN,t] = GSP_MakeGLM({X},[ones(395,1),Regressor'],[0;1]);
+
+% Indices of the regions to retain
+idx_GLM = find(SIGN == 1);
+
+% Bandwidth to use for the Slepian design; for now we go for something
+% large enough
+BW = 838;
+
+% Number of regions
+n_ROI = size(W,1);
+
+% Computation of Slepian vectors
+[Slep,Slep_coefs,Xi,Mu] = GSP_ComputeSlepians(W,BW,idx_GLM,{X});
+
+PlotBrainGraph(zeros(n_ROI,n_ROI),abs(Slep(:,sidx)),Slep(:,sidx),CodeBook,0.5,max(abs(Slep(:,sidx))),...
+ max(abs(Slep(:,sidx))),2,1,'jet','jet',1,...
+ 0.6,[]);
+
+X_SLEP{1} = Slep' * X;
+
+% Computation of Slepian vectors using a random set of nodes as a selection
+% idx_random = randperm(n_ROI);
+% [Slep_rand,Slep_coefs_rand,Xi_rand,Mu_rand] = GSP_ComputeSlepians(W,BW,idx_random(1:length(idx_GLM)),{X});
+
+TB1 = figure;
+set(gcf,'color','w');
+imagesc(X_SLEP{1});
+caxis([-100 100]);
+
+% TB2 = figure;
+% set(gcf,'color','w');
+% imagesc(X_slep_rand{1});
+% caxis([-100 100]);
+
+% Options for the decomposition: we have a real symmetric matrix at
+% hand
+opts.issym=1;
+opts.isreal=1;
+opts.maxit=2500;
+opts.disp=0;
+
+% U contains the eigenvectors, while Lambda is a diagonal matrix
+% containing the eigenvalues of the decomposition. They are sorted in
+% increasing order of amplitude
+[V,Lambda] = eigs(diag(sum(W,2))-W,n_ROI,'SA',opts);
+
+[V_loc,Lambda_loc] = eigs(diag(sum(W(idx_GLM,idx_GLM),2))-W(idx_GLM,idx_GLM),length(idx_GLM),'SA',opts);
+
+X_GFT{1} = V' * X;
+
+X_GFT_loc{1} = V_loc' * X(idx_GLM,:);
+
+TB3 = figure;
+set(gcf,'color','w');
+imagesc(X_GFT{1});
+caxis([-100 100]);
+
+% Number of null realisations to create
+n_null = 100;
+
+% Number of subjects (for now, 1)
+n_subjects = 1;
+
+
+% Generation of my null data
+for s = 1:n_subjects
+
+ disp(num2str(s));
+
+ for n = 1:n_null
+
+ % Filter that we will use for null data generation
+ Phi_2 = zeros(n_ROI,n_ROI);
+ for i = 1:n_ROI
+ tmp = randn;
+
+ if tmp > 0
+ Phi_2(i,i) = 1;
+ else
+ Phi_2(i,i) = -1;
+ end
+ end
+
+ % Filter that we will use for null data generation (Slepians)
+ Phi_3 = zeros(BW,BW);
+ for i = 1:BW
+ tmp = randn;
+
+ if tmp > 0
+ Phi_3(i,i) = 1;
+ else
+ Phi_3(i,i) = -1;
+ end
+ end
+
+ XNull_GFT{s}(:,:,n) = V*Phi_2*V'*X;
+ XNull_SLEP{s}(:,:,n) = Slep*Phi_3*Slep'*X;
+ end
+
+ TNull_GFT{s} = prctile(XNull_GFT{s},[2.5 97.5],3);
+ TNull_SLEP{s} = prctile(XNull_SLEP{s},[2.5 97.5],3);
+
+ EXC_GFT{s} = zeros(n_ROI,395);
+ EXC_GFT{s}(X > TNull_GFT{1}(:,:,2)) = 1;
+ EXC_GFT{s}(X < TNull_GFT{1}(:,:,1)) = -1;
+
+ EXC_SLEP{s} = zeros(n_ROI,395);
+ EXC_SLEP{s}(X > TNull_SLEP{1}(:,:,2)) = 1;
+ EXC_SLEP{s}(X < TNull_SLEP{1}(:,:,1)) = -1;
+end
+
+
+TB4 = figure;
+set(gcf,'color','w');
+imagesc(EXC_GFT{1});
+
+TB5 = figure;
+set(gcf,'color','w');
+imagesc(EXC_SLEP{1});
+
+Slep = -1 * Slep;
+
+% Cutoff on Xi to retain Slepians
+epsilon = 1e-2;
+
+% Computation of concentration
+for t = 1:395
+
+ C(t) = ((Slep(:,diag(Xi)>epsilon)'*X(:,t)))' * Mu(diag(Xi)>epsilon)';
+
+ for n = 1:n_null
+ CN(t,n) = ((Slep(:,diag(Xi)>epsilon)'*squeeze(XNull_SLEP{1}(:,t,n))))' * Mu(diag(Xi)>epsilon)';
+ end
+end
+
+figure;
+set(gcf,'color','w');
+plot_ci(time',[prctile(CN,2.5,2),prctile(CN,97.5,2)],'PatchColor',...
+ [0.2 0.2 0.2], 'PatchAlpha', 0.9, 'MainLineWidth', 2);
+hold on;
+plot(time,C,'r','LineWidth',2);
+plot(time,600*(Regressor_real-mean(Regressor_real)),'color',[1 0.8 0.2]);
+plot(time,600*(Regressor-mean(Regressor)),'color',[0.8 0.4 0]);
+xlabel('Time [s]');
+ylabel('Concentration within visual system');
+xlim([time(1) time(end)]);
+
+
+% Going for trajectory plotting now
+
+
+% Option 1: non linear dimensonality reduction using graphs
+W1 = construct_knn_graph(X_GFT{1}',20,'euclidean');
+W2 = construct_knn_graph(X_SLEP{1}',20,'euclidean');
+W3 = construct_knn_graph(X_SLEP{1}(diag(Xi)>epsilon,:)',20,'euclidean');
+W4 = construct_knn_graph(X_GFT_loc{1}',20,'euclidean');
+
+[X1,Y1,Z1] = compute_non_linear_dim_reduction(W1);
+[X2,Y2,Z2] = compute_non_linear_dim_reduction(W2);
+[X3,Y3,Z3] = compute_non_linear_dim_reduction(W3);
+[X4,Y4,Z4] = compute_non_linear_dim_reduction(W4);
+
+
+CV = Regressor_real;
+figure;
+scatter3(X1,Y1,Z1,10,CV,'filled');
+colormap(jet);
+caxis([min(Regressor_real),max(Regressor_real)]);
+
+CV = Regressor_real;
+figure;
+scatter3(X3,Y3,Z3,10,CV,'filled');
+colormap(jet);
+caxis([min(Regressor_real),max(Regressor_real)]);
+
+CV = Regressor_real;
+figure;
+scatter3(X2,Y2,Z2,10,CV,'filled');
+colormap(jet);
+caxis([min(Regressor_real),max(Regressor_real)]);
+
+CV = Regressor_real;
+figure;
+scatter3(X4,Y4,Z4,10,CV,'filled');
+colormap(jet);
+caxis([min(Regressor_real),max(Regressor_real)]);
+
+
+% Option 2: PCA (linear reduction of dimensionality)
+
+% Data has size N x D
+Data = X_SLEP{1}(diag(Xi)>epsilon,:)';
+
+% Demeaning
+Mu_data = mean(Data,1);
+Data = Data - repmat(Mu_data,[size(Data,1),1]);
+
+% Computation of covariance matrix
+C_Data = 1/size(Data,1) * (Data' * Data);
+
+% Computation of its eigenvectors
+[V_PCA,Lambda_PCA] = eigs(C_Data,size(C_Data,1),'LA',opts);