diff --git a/toolbox/core/bst_colormaps.m b/toolbox/core/bst_colormaps.m
index 8e776066d..3e45747a6 100644
--- a/toolbox/core/bst_colormaps.m
+++ b/toolbox/core/bst_colormaps.m
@@ -142,6 +142,7 @@
             sColormap.CMap             = cmap_viridis(DEFAULT_CMAP_SIZE);
             sColormap.isAbsoluteValues = 0;
             sColormap.MaxMode          = 'global';
+            sColormap.isRealMin        = 1;
         % Time-frequency maps
         case 'timefreq'
             sColormap.Name             = 'cmap_magma';
diff --git a/toolbox/gui/panel_surface.m b/toolbox/gui/panel_surface.m
index 2ab2cf4a7..864b0d88c 100644
--- a/toolbox/gui/panel_surface.m
+++ b/toolbox/gui/panel_surface.m
@@ -1561,7 +1561,11 @@ function UpdateSurfaceProperties()
                 % === CHECKS ===
                 % If min/max values for this file were not computed yet
                 if isempty(TessInfo(iTess).DataMinMax)
-                    TessInfo(iTess).DataMinMax = bst_memory('GetResultsMaximum', iDS, iResult);
+                    if isequal(GlobalData.DataSet(iDS).Results(iResult).ColormapType, 'time')
+                        TessInfo(iTess).DataMinMax = GlobalData.DataSet(iDS).Results(iResult).Time;
+                    else
+                        TessInfo(iTess).DataMinMax = bst_memory('GetResultsMaximum', iDS, iResult);
+                    end
                 end
                 % Reset Overlay cube
                 TessInfo(iTess).OverlayCube = [];
diff --git a/toolbox/io/import_sources.m b/toolbox/io/import_sources.m
index 8e452fee4..79c5dcfcd 100644
--- a/toolbox/io/import_sources.m
+++ b/toolbox/io/import_sources.m
@@ -176,6 +176,7 @@
 %% ===== READ SOURCE FILES =====
 sMri = [];
 maps = cell(1, length(SourceFiles));
+SPM = cell(1, length(SourceFiles));
 % Loop on each input file
 for iFile = 1:length(SourceFiles)
     % === GET FILE TYPE ===  
@@ -241,27 +242,31 @@
         end
         break;
     end
+    % Load SPM results
+    if isStat
+        % Load SPM.mat
+        SpmMat = load(fullfile(fPath, 'SPM.mat'));
+        SPM{iFile} = SpmMat.SPM;
+        % Add sorted T values in the file if importing from a volume (not keeping the full distribution otherwise)
+        if ~isempty(grid)
+            SPM{iFile}.SortedT = sort(sMriSrc.Cube(:));
+        end
+    end
 end
 % Concatenate in time all the selected files
 map = cat(2, maps{:});
 
 % === STATISTICAL THRESHOLD (SPM) ===
 if isStat
-    % Load SPM.mat
-    SpmMat = load(fullfile(fPath, 'SPM.mat'));
     % New stat/results structure
     ResultsMat = db_template('statmat');
     ResultsMat.tmap       = map;
     ResultsMat.pmap       = [];
     ResultsMat.df         = [];
-    ResultsMat.SPM        = SpmMat.SPM;
+    ResultsMat.SPM        = SPM;
     ResultsMat.Correction = 'no';
     ResultsMat.Type       = 'results';
     FileType = 'presults';
-    % Add sorted T values in the file if importing from a volume (not keeping the full distribution otherwise)
-    if ~isempty(grid)
-        ResultsMat.SPM.SortedT = sort(sMriSrc.Cube(:));
-    end
     % Time vector
     if isempty(TimeVector) || (length(TimeVector) ~= size(ResultsMat.tmap,2))
         ResultsMat.Time = 0:(size(map,2)-1);
diff --git a/toolbox/process/functions/process_epileptogenicity.m b/toolbox/process/functions/process_epileptogenicity.m
index c9e4bcec6..732241d0f 100644
--- a/toolbox/process/functions/process_epileptogenicity.m
+++ b/toolbox/process/functions/process_epileptogenicity.m
@@ -145,6 +145,8 @@
     sMri = in_mri_bst(sSubject.Anatomy(sSubject.iAnatomy).FileName);
 
     % ===== EXPORT INPUT FILES =====
+    % Empty temporary folder
+    gui_brainstorm('EmptyTempFolder');
     % Work in Brainstorm's temporary folder
     workDir = bst_fullfile(bst_get('BrainstormTmpDir'), 'ImaGIN_epileptogenicity');
     % Make sure Matlab is not currently in the work directory
@@ -364,7 +366,7 @@
         % File comment = File name
         [tmp, Comment] = bst_fileparts(listFiles(i).name);
         % Import file
-        tmpFile = import_sources(iStudy, [], bst_fullfile(workDir, listFiles(i).name), [], fileFormat, Comment, 's', [fileLatency(1), fileLatency(end)]);
+        tmpFile = import_sources(iStudy, [], bst_fullfile(workDir, listFiles(i).name), [], fileFormat, Comment, 's', [min(fileLatency), max(fileLatency)]);
         % Add history field
         tmpFile = bst_history('add', tmpFile, 'epilepto', strHistory);
     end
@@ -437,7 +439,3 @@
     end
 end
 
-    
-
-
-
diff --git a/toolbox/process/functions/process_extract_pthresh.m b/toolbox/process/functions/process_extract_pthresh.m
index e028422b9..23ee2fab2 100644
--- a/toolbox/process/functions/process_extract_pthresh.m
+++ b/toolbox/process/functions/process_extract_pthresh.m
@@ -246,7 +246,7 @@
     OutputFiles{1} = DataFile;
 end
 
-
+            
 %% ===== APPLY THRESHOLD =====
 function [threshmap, tThreshUnder, tThreshOver] = Compute(StatMat, StatThreshOptions)
     % If options not provided, read them from the interface
@@ -288,22 +288,18 @@
             warning('SPM must be in the Matlab path to compute the statistical thresold for this file.');
             pmask = ones(size(StatMat.tmap));
         else
-            % Compute threshold for statistical map
-            df = [StatMat.SPM.xCon(1).eidf, StatMat.SPM.xX.erdf];
-            S = StatMat.SPM.xVol.S;    %-search Volume {voxels}
-            R = StatMat.SPM.xVol.R;    %-search Volume {resels}
-            % Correction
-            switch (StatThreshOptions.Correction)
-                % Note: always one-sided t-test, really the case?
-                case {'none', 'no'}
-                    tThreshOver = spm_u(StatThreshOptions.pThreshold, df, 'T');
-                case 'bonferroni'
-                    tThreshOver = spm_uc_Bonf(StatThreshOptions.pThreshold, df, 'T', S, 1);
-                case 'fdr'
-                    tThreshOver = spm_uc(StatThreshOptions.pThreshold, df, 'T', R, 1, S);
+            % Threshold each time point independently
+            if iscell(StatMat.SPM)
+                pmask = false(size(StatMat.tmap));
+                for iTime = 1:size(StatMat.tmap, 2)
+                    tThreshOver = GetSpmThreshold(StatMat.SPM{iTime}, StatThreshOptions.Correction, StatThreshOptions.pThreshold);
+                    pmask(:,iTime) = (StatMat.tmap(:,iTime) >= tThreshOver);
+                end
+            % Threshold all the time points in the same way
+            else
+                tThreshOver = GetSpmThreshold(StatMat.SPM, StatThreshOptions.Correction, StatThreshOptions.pThreshold);
+                pmask = (StatMat.tmap >= tThreshOver);
             end
-            % Activated voxels
-            pmask = (StatMat.tmap >= tThreshOver);
             tThreshUnder = [];
         end
     else
@@ -400,3 +396,25 @@
     iMaxUnder = [];
 end
 end
+
+
+%% ===== COMPUTE SPM THRESHOLD =====
+function tThreshOver = GetSpmThreshold(SPM, Method, pThreshold)
+    % Compute threshold for statistical map
+    df = [SPM.xCon(1).eidf, SPM.xX.erdf];
+    S = SPM.xVol.S;    %-search Volume {voxels}
+    R = SPM.xVol.R;    %-search Volume {resels}
+    % Correction
+    switch (Method)
+        % Note: always one-sided t-test, really the case?
+        case {'none', 'no'}
+            tThreshOver = spm_u(pThreshold, df, 'T');
+        case 'bonferroni'
+            tThreshOver = spm_uc_Bonf(pThreshold, df, 'T', S, 1);
+        case 'fdr'
+            tThreshOver = spm_uc(pThreshold, df, 'T', R, 1, S);
+    end
+end
+
+
+