diff --git a/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/analyzeMgPipeResults.m b/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/analyzeMgPipeResults.m
index 8bfadd3f13..6ccc05087c 100755
--- a/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/analyzeMgPipeResults.m
+++ b/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/analyzeMgPipeResults.m
@@ -112,9 +112,9 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
             % Print the results as a text file
             filename = strrep(fileList{i},'.csv','');
             filename = strrep(filename,'.txt','');
-            writetable(cell2table(Statistics),[statPath filesep filename '_' sampleGroupHeaders{j} '_Statistics'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+            writetable(cell2table(Statistics),[statPath filesep filename '_' sampleGroupHeaders{j} '_Statistics.csv'],'WriteVariableNames',false);
             if size(significantFeatures,2)>1
-                writetable(cell2table(significantFeatures),[statPath filesep filename '_' sampleGroupHeaders{j} '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+                writetable(cell2table(significantFeatures),[statPath filesep filename '_' sampleGroupHeaders{j} '_SignificantFeatures.csv'],'WriteVariableNames',false);
             end
             
             % create plots
@@ -133,9 +133,9 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
         % Print the results as a text file
         filename = strrep(fileList{i},'.csv','');
         filename = strrep(filename,'.txt','');
-        writetable(cell2table(Statistics),[statPath filesep filename '_Statistics'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+        writetable(cell2table(Statistics),[statPath filesep filename '_Statistics.csv'],'WriteVariableNames',false);
         if size(significantFeatures,2)>1
-            writetable(cell2table(significantFeatures),[statPath filesep filename '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+            writetable(cell2table(significantFeatures),[statPath filesep filename '_SignificantFeatures.csv'],'WriteVariableNames',false);
         end
         
         % create violin plots
diff --git a/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/makeViolinPlots.m b/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/makeViolinPlots.m
index cd0094afa6..6c4db5c06a 100755
--- a/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/makeViolinPlots.m
+++ b/src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/makeViolinPlots.m
@@ -54,9 +54,16 @@ function makeViolinPlots(sampleData, sampleInformation, varargin)
 cnt=1;
 delArray=[];
 for i=2:size(sampleData,1)
-    if sum(str2double(sampleData(i,2:end)))<0.0001
-        delArray(cnt,1)=i;
-        cnt=cnt+1;
+    if contains(version,'R202') % for MATLAB 2020a or newer
+        if sum(cell2mat(sampleData(i,2:end)))<0.0001
+            delArray(cnt,1)=i;
+            cnt=cnt+1;
+        end
+    else
+        if sum(str2double(sampleData(i,2:end)))<0.0001
+            delArray(cnt,1)=i;
+            cnt=cnt+1;
+        end
     end
 end
 sampleData(delArray,:)=[];
diff --git a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/buildModelStorage.m b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/buildModelStorage.m
index 8dd0f6d006..ee00f69b77 100755
--- a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/buildModelStorage.m
+++ b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/buildModelStorage.m
@@ -50,6 +50,8 @@
 % define human-derived metabolites present in the gut: primary bile acids, amines, mucins, host glycans
 if includeHumanMets
     HumanMets={'gchola','-10';'tdchola','-10';'tchola','-10';'dgchol','-10';'34dhphe','-10';'5htrp','-10';'Lkynr','-10';'f1a','-1';'gncore1','-1';'gncore2','-1';'dsT_antigen','-1';'sTn_antigen','-1';'core8','-1';'core7','-1';'core5','-1';'core4','-1';'ha','-1';'cspg_a','-1';'cspg_b','-1';'cspg_c','-1';'cspg_d','-1';'cspg_e','-1';'hspg','-1'};
+else
+    HumanMets={};
 end
 
 % get all exchanges that can carry flux in at least one model on the given
@@ -94,11 +96,14 @@
             end
         end
         model = useDiet(model, diet,0);
-
+        
         if includeHumanMets
-            % add the human metabolites
+            % add the human metabolites if not already included
+            % in the diet
             for l=1:length(HumanMets)
-                model=changeRxnBounds(model,strcat('EX_',HumanMets{l},'(e)'),str2num(HumanMets{l,2}),'l');
+                if isempty(find(strcmp(diet(:,1),['EX_',HumanMets{l},'(e)'])))
+                    model=changeRxnBounds(model,['EX_',HumanMets{l},'(e)'],str2num(HumanMets{l,2}),'l');
+                end
             end
         end
     end
@@ -174,11 +179,14 @@
                 end
             end
             model = useDiet(model, diet,0);
-
+            
             if includeHumanMets
-                % add the human metabolites
+                % add the human metabolites if not already included
+                % in the diet
                 for l=1:length(HumanMets)
-                    model=changeRxnBounds(model,strcat('EX_',HumanMets{l},'(e)'),str2num(HumanMets{l,2}),'l');
+                    if isempty(find(strcmp(diet(:,1),['EX_',HumanMets{l},'(e)'])))
+                        model=changeRxnBounds(model,['EX_',HumanMets{l},'(e)'],str2num(HumanMets{l,2}),'l');
+                    end
                 end
             end
         end
diff --git a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/createPanModels.m b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/createPanModels.m
index 6740d11b88..a6810ac78a 100755
--- a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/createPanModels.m
+++ b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/createPanModels.m
@@ -1,4 +1,4 @@
-function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorkers)
+function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorkers, builtTaxa)
 % This function creates pan-models for all unique taxa (e.g., species)
 % included in the AGORA resource. If reconstructions of multiple strains
 % in a given taxon are present, the reactions in these reconstructions will
@@ -23,16 +23,22 @@ function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorker
 %                  created pan-models will be stored in. Must end with a file separator.
 %    taxonLevel    String with desired taxonomical level of the pan-models.
 %                  Allowed inputs are 'Species','Genus','Family','Order', 'Class','Phylum'.
-%    agoraVersion  Version of AGORA that will be used (allowed inputs: 'AGORA', 'AGORA2')
+%    agoraVersion  Version of AGORA that will be used (allowed inputs: 'AGORA', 'AGORA2', 
+%                  alternatively: path to custom table with reconstruction information)
 %
 % OPTIONAL INPUTS
 %    numWorkers    Number of workers for parallel pool (default: no pool)
+%    builtTaxa     Names of taxa in table that will be built (default:
+%                  all). Need to be entered as a cell array of strings with names written
+%                  exactly as in the corresponding column in the table.
 %
 % .. Authors
 %       - Stefania Magnusdottir, 2016
 %       - Almut Heinken, 06/2018: adapted to function.
 %       - Almut Heinken, 03/2021: enabled parallelization
 %       - Almut Heinken, 05/2024: enforced specifying AGORA version
+%       - Almut Heinken, 10/2024: allowed specifying input table and taxa
+%                                 to create
 
 tol = 1e-5;
 
@@ -58,13 +64,13 @@ function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorker
 solver = CBT_LP_SOLVER;
 environment = getEnvironment();
 
-% load the file with information on AGORA/AGORA2
+% load the file with information on AGORA/AGORA2 or custom table
 if strcmp(agoraVersion,'AGORA')
     infoFile = readInputTableForPipeline('AGORA_infoFile.xlsx');
 elseif strcmp(agoraVersion,'AGORA2')
     infoFile = readInputTableForPipeline('AGORA2_infoFile.xlsx');
 else
-    error('Incorrect input for agoraVersion! Allowed inputs: AGORA, AGORA2')
+    infoFile = readInputTableForPipeline(agoraVersion);
 end
 
 % get the reaction and metabolite database
@@ -78,6 +84,12 @@ function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorker
 allTaxa(strncmp(allTaxa, 'unclassified',  12)) = [];
 allTaxa(strcmp(allTaxa, '')) = [];
 
+% reduce to specified list of taxa to build if entered
+if nargin >5
+    [~,I] = setdiff(allTaxa,builtTaxa);
+    allTaxa(I,:) = [];
+end
+
 % Remove models that have already been assembled from the list of models to create
 dInfo = dir(panPath);
 dInfo = dInfo(~[dInfo.isdir]);
diff --git a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/microbiotaModelSimulator.m b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/microbiotaModelSimulator.m
old mode 100644
new mode 100755
index 405aacf285..11f8a30b50
--- a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/microbiotaModelSimulator.m
+++ b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/microbiotaModelSimulator.m
@@ -79,7 +79,11 @@
 allDietExch = regexprep(allDietExch,'\[fe\]','\[d\]');
 
 % define human-derived metabolites present in the gut: primary bile acids, amines, mucins, host glycans
-HumanMets={'gchola','-10';'tdchola','-10';'tchola','-10';'dgchol','-10';'34dhphe','-10';'5htrp','-10';'Lkynr','-10';'f1a','-1';'gncore1','-1';'gncore2','-1';'dsT_antigen','-1';'sTn_antigen','-1';'core8','-1';'core7','-1';'core5','-1';'core4','-1';'ha','-1';'cspg_a','-1';'cspg_b','-1';'cspg_c','-1';'cspg_d','-1';'cspg_e','-1';'hspg','-1'};
+if includeHumanMets
+    HumanMets={'gchola','-10';'tdchola','-10';'tchola','-10';'dgchol','-10';'34dhphe','-10';'5htrp','-10';'Lkynr','-10';'f1a','-1';'gncore1','-1';'gncore2','-1';'dsT_antigen','-1';'sTn_antigen','-1';'core8','-1';'core7','-1';'core5','-1';'core4','-1';'ha','-1';'cspg_a','-1';'cspg_b','-1';'cspg_c','-1';'cspg_d','-1';'cspg_e','-1';'hspg','-1'};
+else
+    HumanMets = {};
+end
 
 %% start the simulations
 
@@ -267,7 +271,7 @@
                     model=changeRxnBounds(model,['Host_' hostBiomassRxn],0.001,'l');
                     model=changeRxnBounds(model,['Host_' hostBiomassRxn],hostBiomassRxnFlux,'u');
                 end
-
+                
                 solution_allOpen = optimizeCbModel(model);
                 % solution_allOpen=solveCobraLPCPLEX(model,2,0,0,[],0);
                 if solution_allOpen.stat==0
@@ -281,13 +285,13 @@
                     FecalRxn = AllRxn(FecalInd);
                     FecalRxn=setdiff(FecalRxn,'EX_microbeBiomass[fe]','stable');
                     DietRxn = AllRxn(DietInd);
-
+                    
                     %% computing fluxes on the rich diet
                     if rDiet==1 && computeProfiles
                         % remove exchanges that cannot carry flux
                         FecalRxn=intersect(FecalRxn,allFecalExch);
                         DietRxn=intersect(DietRxn,allDietExch);
-
+                        
                         [minFlux,maxFlux]=guidedSim(model,FecalRxn);
                         minFluxFecal = minFlux;
                         maxFluxFecal = maxFlux;
@@ -304,32 +308,35 @@
                             netUptakeTmp{k}{1}{index,3} = minFluxFecal(i,1);
                         end
                     end
-
+                    
                     %% Computing fluxes on the input diet
-
+                    
                     % remove exchanges that cannot carry flux
                     FecalRxn=intersect(FecalRxn,allFecalExch);
                     DietRxn=intersect(DietRxn,allDietExch);
-
+                    
                     model_sd=model;
                     if adaptMedium
                         [diet] = adaptVMHDietToAGORA(loadDiet,'Microbiota');
                     else
                         diet = readInputTableForPipeline(loadDiet);  % load the text file with the diet
- 
-                        for j = 2:length(diet)
+                        
+                        for j = 1:length(diet)
                             diet{j, 2} = num2str(-(diet{j, 2}));
                         end
                     end
                     [model_sd] = useDiet(model_sd, diet,0);
-
+                    
                     if includeHumanMets
-                        % add the human metabolites
+                        % add the human metabolites if not already included
+                        % in the diet
                         for l=1:length(HumanMets)
-                            model_sd=changeRxnBounds(model_sd,strcat('Diet_EX_',HumanMets{l},'[d]'),str2num(HumanMets{l,2}),'l');
+                            if isempty(find(strcmp(diet(:,1),['Diet_EX_',HumanMets{l},'[d]'])))
+                                model_sd=changeRxnBounds(model_sd,['Diet_EX_',HumanMets{l},'[d]'],str2num(HumanMets{l,2}),'l');
+                            end
                         end
                     end
-
+                    
                     solution_sDiet=optimizeCbModel(model_sd);
                     % solution_sDiet=solveCobraLPCPLEX(model_sd,2,0,0,[],0);
                     growthRatesTmp{k}{2}=solution_sDiet.f;
@@ -356,64 +363,64 @@
                                 netUptakeTmp{k}{2}{index,3} = minFluxFecal(i,1);
                             end
                         end
-
-                        microbiota_model=model_sd;
-                        parsave([resPath filesep 'Diet' filesep 'microbiota_model_diet_' sampleID '.mat'],microbiota_model)
-
-                        %% Using personalized diet not documented in MgPipe and bug checked yet!!!!
-
-                        if pDiet==1
-                            model_pd=model;
-                            [Numbers, Strings] = xlsread(strcat(abundancepath,fileNameDiets));
-                            % diet exchange reactions
-                            DietNames = Strings(2:end,1);
-                            % Diet exchanges for all individuals
-                            Diets(:,k) = cellstr(num2str((Numbers(1:end,k))));
-                            Dietexchanges = {DietNames{:,1} ; Diets{:,k}}';
-                            Dietexchanges = regexprep(Dietexchanges,'EX_','Diet_EX_');
-                            Dietexchanges = regexprep(Dietexchanges,'\(e\)','\[d\]');
-
-                            model_pd = setDietConstraints(model_pd,Dietexchanges);
-
-                            if includeHumanMets
-                                % add the human metabolites
-                                for l=1:length(HumanMets)
-                                    model_pd=changeRxnBounds(model_pd,strcat('Diet_EX_',HumanMets{l},'[d]'),str2num(HumanMets{l,2}),'l');
-                                end
+                    end
+                    
+                    microbiota_model=model_sd;
+                    parsave([resPath filesep 'Diet' filesep 'microbiota_model_diet_' sampleID '.mat'],microbiota_model)
+                    
+                    %% Using personalized diet not documented in MgPipe and bug checked yet!!!!
+                    
+                    if pDiet==1
+                        model_pd=model;
+                        [Numbers, Strings] = xlsread(strcat(abundancepath,fileNameDiets));
+                        % diet exchange reactions
+                        DietNames = Strings(2:end,1);
+                        % Diet exchanges for all individuals
+                        Diets(:,k) = cellstr(num2str((Numbers(1:end,k))));
+                        Dietexchanges = {DietNames{:,1} ; Diets{:,k}}';
+                        Dietexchanges = regexprep(Dietexchanges,'EX_','Diet_EX_');
+                        Dietexchanges = regexprep(Dietexchanges,'\(e\)','\[d\]');
+                        
+                        model_pd = setDietConstraints(model_pd,Dietexchanges);
+                        
+                        if includeHumanMets
+                            % add the human metabolites
+                            for l=1:length(HumanMets)
+                                model_pd=changeRxnBounds(model_pd,strcat('Diet_EX_',HumanMets{l},'[d]'),str2num(HumanMets{l,2}),'l');
                             end
-
-                            solution_pdiet=optimizeCbModel(model_pd);
-                            %solution_pdiet=solveCobraLPCPLEX(model_pd,2,0,0,[],0);
-                            growthRatesTmp{k}{3}=solution_pdiet.f;
-                            if solution_pdiet.stat==0
-                                warning('growthRates detected one or more infeasible models. Please check infeasModels object !')
-                                infeasModelsTmp{k} = model.name;
-                                netProductionTmp{k}{3} = {};
-                                netUptakeTmp{k}{3} = {};
-                            else
-                                if computeProfiles
-                                    [minFlux,maxFlux]=guidedSim(model_pd,FecalRxn);
-                                    minFluxFecal = minFlux;
-                                    maxFluxFecal = maxFlux;
-                                    [minFlux,maxFlux]=guidedSim(model_pd,DietRxn);
-                                    minFluxDiet = minFlux;
-                                    maxFluxDiet = maxFlux;
-                                    netProductionTmp{k}{3}=exchanges;
-                                    netUptakeTmp{k}{3}=exchanges;
-                                    for i =1:length(FecalRxn)
-                                        [truefalse, index] = ismember(FecalRxn(i), exchanges);
-                                        netProductionTmp{k}{3}{index,2} = minFluxDiet(i,1);
-                                        netProductionTmp{k}{3}{index,3} = maxFluxFecal(i,1);
-                                        netUptakeTmp{k}{3}{index,2} = maxFluxDiet(i,1);
-                                        netUptakeTmp{k}{3}{index,3} = minFluxFecal(i,1);
-                                    end
+                        end
+                        
+                        solution_pdiet=optimizeCbModel(model_pd);
+                        %solution_pdiet=solveCobraLPCPLEX(model_pd,2,0,0,[],0);
+                        growthRatesTmp{k}{3}=solution_pdiet.f;
+                        if solution_pdiet.stat==0
+                            warning('growthRates detected one or more infeasible models. Please check infeasModels object !')
+                            infeasModelsTmp{k} = model.name;
+                            netProductionTmp{k}{3} = {};
+                            netUptakeTmp{k}{3} = {};
+                        else
+                            if computeProfiles
+                                [minFlux,maxFlux]=guidedSim(model_pd,FecalRxn);
+                                minFluxFecal = minFlux;
+                                maxFluxFecal = maxFlux;
+                                [minFlux,maxFlux]=guidedSim(model_pd,DietRxn);
+                                minFluxDiet = minFlux;
+                                maxFluxDiet = maxFlux;
+                                netProductionTmp{k}{3}=exchanges;
+                                netUptakeTmp{k}{3}=exchanges;
+                                for i =1:length(FecalRxn)
+                                    [truefalse, index] = ismember(FecalRxn(i), exchanges);
+                                    netProductionTmp{k}{3}{index,2} = minFluxDiet(i,1);
+                                    netProductionTmp{k}{3}{index,3} = maxFluxFecal(i,1);
+                                    netUptakeTmp{k}{3}{index,2} = maxFluxDiet(i,1);
+                                    netUptakeTmp{k}{3}{index,3} = minFluxFecal(i,1);
                                 end
-
-                                % save the model with personalized diet
-                                microbiota_model=model_pd;
-                                mkdir(strcat(resPath,'Personalized'))
-                                parsave([resPath filesep 'Personalized' filesep 'microbiota_model_pDiet_' sampleID '.mat'],microbiota_model)
                             end
+                            
+                            % save the model with personalized diet
+                            microbiota_model=model_pd;
+                            mkdir(strcat(resPath,'Personalized'))
+                            parsave([resPath filesep 'Personalized' filesep 'microbiota_model_pDiet_' sampleID '.mat'],microbiota_model)
                         end
                     end
                 end
@@ -435,14 +442,14 @@
                 end
             end
             if ~isempty(growthRatesTmp{k})
-                 if ~isempty(growthRatesTmp{k}{1})
+                if ~isempty(growthRatesTmp{k}{1})
                     growthRates{k+1,2} = growthRatesTmp{k}{1};
                     growthRates{k+1,3} = growthRatesTmp{k}{2};
                     if length(growthRatesTmp{k})>2
                         growthRates{1,4} = 'Personalized diet';
                         growthRates{k+1,4} = growthRatesTmp{k}{3};
                     end
-                 end
+                end
             end
             if ~isempty(infeasModelsTmp) && k <= length(infeasModelsTmp)
                 infeasModels{k,1} = infeasModelsTmp{k};
diff --git a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/retrieveModelStats.m b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/retrieveModelStats.m
index eecc1bc0d4..dea6f321df 100644
--- a/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/retrieveModelStats.m
+++ b/src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/retrieveModelStats.m
@@ -101,15 +101,18 @@
 % create violin plot of model stats
 if nargin <5
     % have reactions and metabolites in one plot
-    figure
-    subplot(1,2,1)
-    violinplot(data(:,1:2),{'Reactions','Metabolites'});
-    set(gca, 'FontSize', 12)
-    subplot(1,2,2)
-    violinplot(data(:,3),{'Microbes'});
-    set(gca, 'FontSize', 12)
-    sgtitle('Reaction, metabolite and microbe numbers in microbiome models')
-    print('MicrobiomeModel_Sizes','-dpng','-r300')
+    % does not work if all data points are the same
+    if ~numel(unique(data(:,1)))==1 && ~numel(unique(data(:,2)))==1 && ~numel(unique(data(:,3)))==1
+        figure
+        subplot(1,2,1)
+        violinplot(data(:,1:2),{'Reactions','Metabolites'});
+        set(gca, 'FontSize', 12)
+        subplot(1,2,2)
+        violinplot(data(:,3),{'Microbes'});
+        set(gca, 'FontSize', 12)
+        sgtitle('Reaction, metabolite and microbe numbers in microbiome models')
+        print('MicrobiomeModel_Sizes','-dpng','-r300')
+    end
     
 else
     % perform statistical analysis if file with stratification is provided
@@ -180,9 +183,12 @@
     violinplot(data(:,2),infoFile(:,2));
     title('Metabolites')
     set(gca, 'FontSize', 12)
-    subplot(1,3,3)
-    violinplot(data(:,3),infoFile(:,2));
-    title('Microbes')
+    % does not work if all data points are the same
+    if ~numel(unique(data(:,3)))==1
+        subplot(1,3,3)
+        violinplot(data(:,3),infoFile(:,2));
+        title('Microbes')
+    end
     set(gca, 'FontSize', 12)
     hold on
     sgtitle('Reaction, metabolite and microbe numbers in microbiome models')
diff --git a/src/reconstruction/demeter/src/debugging/runDebuggingTools.m b/src/reconstruction/demeter/src/debugging/runDebuggingTools.m
index e08d7326c3..7b522aa5bb 100755
--- a/src/reconstruction/demeter/src/debugging/runDebuggingTools.m
+++ b/src/reconstruction/demeter/src/debugging/runDebuggingTools.m
@@ -94,7 +94,8 @@
     failedModels = union(failedModels,tooHighATP);
 end
 if isfile([testResultsFolder filesep reconVersion '_refined' filesep 'ATP_from_O2_' reconVersion '.txt'])
-    atpFromO2 = table2cell(readtable([testResultsFolder filesep reconVersion '_refined' filesep 'ATP_from_O2_' reconVersion '.txt']));
+    atpFromO2 = readInputTableForPipeline([testResultsFolder filesep reconVersion '_refined' filesep 'ATP_from_O2_' reconVersion '.txt']);
+    atpFromO2{1,2} = str2double(atpFromO2{1,2});
     atpModels = atpFromO2(find(cell2mat(atpFromO2(:,2)) > 0.00001),1);
     failedModels = union(failedModels,atpModels);
 end
diff --git a/src/reconstruction/demeter/src/refinement/createPeriplasmaticSpace.m b/src/reconstruction/demeter/src/refinement/createPeriplasmaticSpace.m
old mode 100644
new mode 100755
index c30b7666b6..4f40f3dd53
--- a/src/reconstruction/demeter/src/refinement/createPeriplasmaticSpace.m
+++ b/src/reconstruction/demeter/src/refinement/createPeriplasmaticSpace.m
@@ -71,6 +71,9 @@
                 transpRxns = findRxnsFromMets(model,exMets{i});
                 % remove exchange reactions
                 transpRxns(find(strncmp(transpRxns,'EX_',3)))=[];
+                % remove reactions that are already periplasmatic
+                forms = printRxnFormula(model,transpRxns);
+                transpRxns(find(contains(forms,'[p]')))=[];
                 for j=1:length(transpRxns)
                     rxnsToAdd{end+1} = [transpRxns{j} 'pp'];
                     rxnNames{end+1} = [model.rxnNames{find(strcmp(model.rxns,transpRxns{j}))} ', periplasmatic'];
diff --git a/src/reconstruction/demeter/src/refinement/performDataDrivenRefinement.m b/src/reconstruction/demeter/src/refinement/performDataDrivenRefinement.m
index 678a4ab4c1..029ebd7b28 100644
--- a/src/reconstruction/demeter/src/refinement/performDataDrivenRefinement.m
+++ b/src/reconstruction/demeter/src/refinement/performDataDrivenRefinement.m
@@ -61,29 +61,26 @@
 %% gapfill if there are any false negatives
 osenseStr='min';
 
-dataDrivenGapfill={};
 if ~isempty(FNs)
     for j=1:length(FNs)
         metExch=['EX_' database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j})),1} '(e)'];
         % find reactions that could be gap-filled to enable flux
         % seems to try to fix exchanges that are not part of the model,
         % leading it to crash:
-        %         Objective reactions not found in model!
+        %  Objective reactions not found in model!
         %
         % Error in runGapfillingFunctions (line 41)
         % model = changeObjective(model, objectiveFunction);
-     %   add exchange reaction
-     if isempty(find(strcmp(model.rxns,metExch)))
-         met = [database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j}))} '(e)'];
-         model = addExchangeRxn(model,met,-1,1000);
-     end
-     [model,gapfilledRxns] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
-     dataDrivenGapfill=union(dataDrivenGapfill,gapfilledRxns);
-     %  end
-    end
-    if ~isempty(dataDrivenGapfill)
-        summary.('DataDrivenGapfill')=dataDrivenGapfill;
+        %   add exchange reaction
+        if isempty(find(strcmp(model.rxns,metExch)))
+            met = [database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j}))} '[e]'];
+            model = addExchangeRxn(model,met,-1000,1000);
+        end
+        [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
     end
+    summary.('conditionSpecificGapfill') = union(summary.('conditionSpecificGapfill'),condGF);
+    summary.('targetedGapfill') = union(summary.('targetedGapfill'),targetGF);
+    summary.('relaxFBAGapfill') = union(summary.('relaxFBAGapfill'),relaxGF);
 end
 
 % Fermentation products
@@ -113,10 +110,11 @@
             metExch = FNs{j};
         end
         % find reactions that could be gap-filled to enable flux
-        try
-        [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
-        catch
+        if isempty(find(strcmp(model.rxns,metExch)))
+            met = [database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j}))} '[e]'];
+            model = addExchangeRxn(model,met,-1000,1000);
         end
+        [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
         summary.('conditionSpecificGapfill') = union(summary.('conditionSpecificGapfill'),condGF);
         summary.('targetedGapfill') = union(summary.('targetedGapfill'),targetGF);
         summary.('relaxFBAGapfill') = union(summary.('relaxFBAGapfill'),relaxGF);
diff --git a/src/reconstruction/demeter/src/refinement/runDemeter.m b/src/reconstruction/demeter/src/refinement/runDemeter.m
index 70b7a0f989..f6168f0001 100644
--- a/src/reconstruction/demeter/src/refinement/runDemeter.m
+++ b/src/reconstruction/demeter/src/refinement/runDemeter.m
@@ -1,4 +1,4 @@
-function [reconVersion,refinedFolder,translatedDraftsFolder,summaryFolder] = runDEMETER(draftFolder, varargin)
+function [reconVersion,refinedFolder,translatedDraftsFolder,summaryFolder] = runDemeter(draftFolder, varargin)
 % This function runs the DEMETER pipeline consisting of three steps:
 % 1) refining all draft reconstructions, 2) testing the refined
 % reconstructions against the input data, 3) preparing a report
@@ -179,9 +179,8 @@
     modelsTmp = {};
     draftModelsTmp = {};
     summariesTmp = {};
-   
-  parfor j=i:i+endPnt
- % for j=i:i+endPnt
+    
+    parfor j=i:i+endPnt
         restoreEnvironment(environment);
         changeCobraSolver(solver, 'LP', 0, -1);
         
diff --git a/test/verifiedTests/analysis/testMultiSpeciesModelling/testCreatePanModels.m b/test/verifiedTests/analysis/testMultiSpeciesModelling/testCreatePanModels.m
index dfa0f901cb..acfb26559a 100644
--- a/test/verifiedTests/analysis/testMultiSpeciesModelling/testCreatePanModels.m
+++ b/test/verifiedTests/analysis/testMultiSpeciesModelling/testCreatePanModels.m
@@ -14,19 +14,33 @@
 fileDir = fileparts(which('testCreatePanModels'));
 cd(fileDir);
 
-% load the AGORA models
-websave('AGORA-master.zip','https://github.com/VirtualMetabolicHuman/AGORA/archive/master.zip')
-try
-    unzip('AGORA-master')
-end
-modPath = [pwd filesep 'AGORA-master' filesep 'CurrentVersion' filesep 'AGORA_1_03' filesep' 'AGORA_1_03_mat'];
-
+modelList={
+        'Abiotrophia_defectiva_ATCC_49176'
+        'Acidaminococcus_fermentans_DSM_20731'
+        'Acidaminococcus_intestini_RyC_MR95'
+        'Acidaminococcus_sp_D21'
+        'Acinetobacter_baumannii_AB0057'
+        'Acinetobacter_calcoaceticus_PHEA_2'
+        'Acinetobacter_haemolyticus_NIPH_261'
+        'Acinetobacter_johnsonii_SH046'
+        'Acinetobacter_junii_SH205'
+        'Acinetobacter_lwoffii_WJ10621'
+        'Acinetobacter_pittii_ANC_4052'
+        'Acinetobacter_radioresistens_NIPH_2130'
+        };
+    for i=1:length(modelList)
+        model = getDistributedModel([modelList{i} '.mat']);
+        % Save all the models into the modelFolder
+        save(fullfile(fileDir, modelList{i}), 'model');
+    end
+    
 numWorkers=4;
 
 % create the pan-models on species level
 panPath=[pwd filesep 'panSpeciesModels'];
 
-createPanModels(modPath,panPath,'Species','AGORA',numWorkers);
+builtTaxa = {'Abiotrophia defectiva','Acidaminococcus fermentans','Acidaminococcus intestini','Acidaminococcus sp. D21','Acinetobacter_baumannii','Acinetobacter haemolyticus','Acinetobacter johnsonii','Acinetobacter junii','Acinetobacter lwoffii','Acinetobacter pittii','Acinetobacter radioresistens'};
+createPanModels(fileDir,panPath,'Species','AGORA_infoFile.xlsx',numWorkers,builtTaxa)
 
 % test that pan-models can grow
 [notGrowing,Biomass_fluxes] = plotBiomassTestResults(panPath, 'pan-models','numWorkers',numWorkers);
@@ -40,7 +54,8 @@
 % create the pan-models on genus level
 panPath=[pwd filesep 'panGenusModels'];
 
-createPanModels(modPath,panPath,'Genus','AGORA',numWorkers);
+builtTaxa = {'Abiotrophia','Acidaminococcus','Acinetobacter'};
+createPanModels(fileDir,panPath,'Genus','AGORA_infoFile.xlsx',numWorkers,builtTaxa);
 
 % test that pan-models can grow
 [notGrowing,Biomass_fluxes] = plotBiomassTestResults(panPath, 'pan-models','numWorkers',numWorkers);
diff --git a/test/verifiedTests/analysis/testMultiSpeciesModelling/testJoinModelsPairwiseFromList.m b/test/verifiedTests/analysis/testMultiSpeciesModelling/testJoinModelsPairwiseFromList.m
index fa3d758e36..758b456623 100644
--- a/test/verifiedTests/analysis/testMultiSpeciesModelling/testJoinModelsPairwiseFromList.m
+++ b/test/verifiedTests/analysis/testMultiSpeciesModelling/testJoinModelsPairwiseFromList.m
@@ -28,8 +28,8 @@
         'Abiotrophia_defectiva_ATCC_49176'
         'Acidaminococcus_fermentans_DSM_20731'
         'Acidaminococcus_intestini_RyC_MR95'
-        'Acidaminococcus_sp_D21'
         'Acinetobacter_calcoaceticus_PHEA_2'
+        'Acinetobacter_baumannii_AB0057'
         };
     for i=1:length(modelList)
         model = getDistributedModel([modelList{i} '.mat']);
diff --git a/test/verifiedTests/analysis/testMultiSpeciesModelling/testSimulatePairwiseInteractions.m b/test/verifiedTests/analysis/testMultiSpeciesModelling/testSimulatePairwiseInteractions.m
index 4d3e4f0b0a..0e2d4fb36c 100644
--- a/test/verifiedTests/analysis/testMultiSpeciesModelling/testSimulatePairwiseInteractions.m
+++ b/test/verifiedTests/analysis/testMultiSpeciesModelling/testSimulatePairwiseInteractions.m
@@ -30,8 +30,8 @@
         'Abiotrophia_defectiva_ATCC_49176'
         'Acidaminococcus_fermentans_DSM_20731'
         'Acidaminococcus_intestini_RyC_MR95'
-        'Acidaminococcus_sp_D21'
         'Acinetobacter_calcoaceticus_PHEA_2'
+        'Acinetobacter_baumannii_AB0057'
         };
     for i=1:length(modelList)
         model = getDistributedModel([modelList{i} '.mat']);