diff --git a/core/FSEOF.m b/core/FSEOF.m
index c07b468d..094766ec 100755
--- a/core/FSEOF.m
+++ b/core/FSEOF.m
@@ -23,7 +23,7 @@
 %
 %   Usage: targets=FSEOF(model,biomassRxn,targetRxn,iterations,coefficient,outputFile)
 %
-%   Eduard Kerkhoven, 2017-12-21
+%   Simonas Marcisauskas, 2018-03-18
 
 if nargin<4
     iterations=10;
@@ -98,9 +98,9 @@
 formatSpec='%s\t%s\t%s\t%s\t%s\t%s\t%s\n';
 if output == 1    %Output to a file
     fid=fopen(outputFile,'w');
-    fprintf(fid,formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystem','Direction','Gr Rule');
+    fprintf(fid,formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystems','Direction','Gr Rule');
 else              %Output to screen
-    fprintf(formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystem','Direction','Gr Rule');
+    fprintf(formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystems','Direction','Gr Rule');
 end
 
 for num=1:length(fseof.target)
@@ -109,7 +109,7 @@
         A1=num2str(num);                                  %row ID
         A2=char(model.rxns(num));                         %enzyme ID
         A3=char(model.rxnNames(num));                     %enzyme Name
-        A4=char(model.subSystems(num));                   %Subsystem
+        A4=char(strjoin(model.subSystems{num,1},';'));                   %Subsystems
         A5=num2str(model.rev(num)*rxnDirection(num,1));   %reaction Dirction
         A6=char(model.grRules(num));                      %Gr Rule
         if output == 1    %Output to a file
diff --git a/core/checkModelStruct.m b/core/checkModelStruct.m
index 40d19ff2..8f8775c6 100755
--- a/core/checkModelStruct.m
+++ b/core/checkModelStruct.m
@@ -14,7 +14,7 @@ function checkModelStruct(model,throwErrors,trimWarnings)
 %
 %   Usage: checkModelStruct(model,throwErrors,trimWarnings)
 %
-%   Simonas Marcisauskas, 2018-03-17
+%   Simonas Marcisauskas, 2018-03-18
 %
 
 if nargin<2
@@ -149,9 +149,11 @@ function checkModelStruct(model,throwErrors,trimWarnings)
     end
 end
 if isfield(model,'subSystems')
-    if ~iscellstr(model.subSystems)
-        EM='The "subSystems" field must be a cell array';
-        dispEM(EM,throwErrors);
+    for i=1:numel(model.subSystems)
+        if ~iscell(model.subSystems{i,1})
+            EM='The "subSystems" field must be a cell array';
+            dispEM(EM,throwErrors);
+        end
     end
 end
 if isfield(model,'eccodes')
diff --git a/core/sortModel.m b/core/sortModel.m
index 80db9c54..757dcc90 100755
--- a/core/sortModel.m
+++ b/core/sortModel.m
@@ -19,7 +19,7 @@
 %
 %   Usage: model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)
 %
-%   Rasmus Agren, 2014-01-09
+%   Simonas Marcisauskas, 2018-03-19
 %
 
 if nargin<2
@@ -65,10 +65,20 @@
        dispEM(EM);
    end
 
-   subsystems=unique(model.subSystems);
+   subsystemsUnique='';
+   subsystemsConcatenated='';
+   for i=1:numel(model.subSystems)
+        subsystemsConcatenated{i,1}=strjoin(model.subSystems{i,1},';');
+        if ~isempty(model.subSystems{i,1})
+            for j=1:numel(model.subSystems{i,1})
+                subsystemsUnique{numel(subsystemsUnique)+1,1}=model.subSystems{i,1}{1,j};
+            end
+        end
+   end
+   subsystemsUnique=unique(subsystemsUnique);
    for i=1:numel(subsystems)
        %Get all reactions for that subsystem
-       rxns=find(ismember(model.subSystems,subsystems(i)));
+       rxns=find(~cellfun(@isempty,regexp(subsystemsConcatenated,subsystemsUnique(i))));
 
        %Temporarily ignore large subsystems because of inefficient
        %implementation
diff --git a/external/hmr/HMRdatabase2_00.mat b/external/hmr/HMRdatabase2_00.mat
index dc516b66..74f64aa8 100644
Binary files a/external/hmr/HMRdatabase2_00.mat and b/external/hmr/HMRdatabase2_00.mat differ
diff --git a/external/kegg/getRxnsFromKEGG.m b/external/kegg/getRxnsFromKEGG.m
index db971d3f..26344ff7 100755
--- a/external/kegg/getRxnsFromKEGG.m
+++ b/external/kegg/getRxnsFromKEGG.m
@@ -50,7 +50,7 @@
 %
 %   Usage: model=getRxnsFromKEGG(keggPath,keepUndefinedStoich,keepIncomplete,keepGeneral)
 %
-%   Simonas Marcisauskas, 2017-05-02
+%   Simonas Marcisauskas, 2018-03-19
 %
 
 %NOTE: This is how one entry looks in the file
@@ -236,17 +236,18 @@
                     pathway=true;
                   end
     
-                  %Don't do this if the pathway is rn01100 (Metabolic pathways)
-                  if ~strcmp('rn01100',tempStruct.value{addToIndex,1})
+                  %Don't save global or overview pathways. The
+                  %ids for such pathways begin with rn011 or rn012;
+                  if ~strcmp('rn011',tempStruct.value{addToIndex,1}(1:5)) && ~strcmp('rn012',tempStruct.value{addToIndex,1}(1:5))
                     model.rxnMiriams{rxnCounter}=tempStruct;
     
-                    %Also save the subSystems entry as being the first path found
-                    if ~any(model.subSystems{rxnCounter})
-                        if strcmp(tline(14:17),'PATH:')
-                            model.subSystems{rxnCounter}=tline(28:end);
-                        else
-                            model.subSystems{rxnCounter}=tline(22:end);
-                        end
+                    %Also save the subSystems names. For the old KEGG
+                    %format, only the first mentioned subsystem is picked.
+                    %Use the newer KEGG format to fetch all the subsystems;
+                    if strcmp(tline(14:17),'PATH:')
+                        model.subSystems{rxnCounter}=tline(28:end);%The old format;
+                    else
+                        model.subSystems{rxnCounter,1}{1,numel(model.subSystems{rxnCounter,1})+1}=tline(22:end);%The new format;
                     end
                   end
               end
diff --git a/external/kegg/keggRxns.mat b/external/kegg/keggRxns.mat
index ce56c680..ecc8610a 100644
Binary files a/external/kegg/keggRxns.mat and b/external/kegg/keggRxns.mat differ
diff --git a/external/metacyc/getRxnsFromMetaCyc.m b/external/metacyc/getRxnsFromMetaCyc.m
index 584e8c6b..a62ce124 100755
--- a/external/metacyc/getRxnsFromMetaCyc.m
+++ b/external/metacyc/getRxnsFromMetaCyc.m
@@ -44,7 +44,7 @@
 %
 %   Usage: model=getRxnsFromMetaCyc(metacycPath,keepTransportRxns,keepUnbalanced,keepUndetermined)
 %
-%   Hao Wang, 2017-11-15
+%   Simonas Marcisauskas, 2018-03-19
 %
 
 %NOTE: This is how one entry looks in the file
@@ -278,11 +278,7 @@
 
 		          [x, y]=ismember(tline(14:end),pwys);
 		          if x
-	          		if isempty(metaCycRxns.subSystems{rxnCounter})
-  	        			metaCycRxns.subSystems{rxnCounter}=pwyNames{y};
-    	      		else
-      	    			metaCycRxns.subSystems{rxnCounter}=strcat(metaCycRxns.subSystems{rxnCounter},';',pwyNames{y});
-        	  		end
+                    metaCycRxns.subSystems{rxnCounter,1}{1,numel(metaCycRxns.subSystems{rxnCounter,1})+1}=pwyNames{y};
 	    	      end
           end
 
diff --git a/external/metacyc/metaCycRxns.mat b/external/metacyc/metaCycRxns.mat
index 0dd0989a..1c1f5182 100644
Binary files a/external/metacyc/metaCycRxns.mat and b/external/metacyc/metaCycRxns.mat differ
diff --git a/io/exportToExcelFormat.m b/io/exportToExcelFormat.m
index 7d82e233..2e16b18c 100755
--- a/io/exportToExcelFormat.m
+++ b/io/exportToExcelFormat.m
@@ -17,7 +17,7 @@ function exportToExcelFormat(model,filename)
 %
 %   Usage: exportToExcelFormat(model,filename)
 %
-%   Simonas Marcisauskas, 2017-09-12
+%   Simonas Marcisauskas, 2018-03-18
 %
 
 [~, A, B]=fileparts(filename);
@@ -164,7 +164,10 @@ function exportToExcelFormat(model,filename)
 rxnSheet=[rxnSheet rxnMiriams];
 
 if isfield(model,'subSystems')
-    rxnSheet=[rxnSheet model.subSystems];
+    for i=1:numel(model.subSystems)
+        subsystems{i,1}=strjoin(model.subSystems{i,1},';');
+    end
+    rxnSheet=[rxnSheet subsystems];
 else
     rxnSheet=[rxnSheet emptyColumn];
 end
diff --git a/io/exportToTabDelimited.m b/io/exportToTabDelimited.m
index 6a021480..852e8647 100755
--- a/io/exportToTabDelimited.m
+++ b/io/exportToTabDelimited.m
@@ -16,7 +16,7 @@ function exportToTabDelimited(model,path)
 %
 %   Usage: exportToTabDelimited(model,path)
 %
-%   Simonas Marcisauskas, 2018-03-17
+%   Simonas Marcisauskas, 2018-03-18
 
 if nargin<2
 	path='./';
@@ -99,7 +99,7 @@ function exportToTabDelimited(model,path)
     end
    
 	if isfield(model,'subSystems')
-        fprintf(rxnFile,[model.subSystems{i} '\t']);
+        fprintf(rxnFile,[strjoin(model.subSystems{i,1},';') '\t']);
     else
         fprintf(rxnFile,'\t');
 	end
diff --git a/io/importExcelModel.m b/io/importExcelModel.m
index bde6ec9c..9f3ad203 100755
--- a/io/importExcelModel.m
+++ b/io/importExcelModel.m
@@ -403,7 +403,7 @@
         case 9
         	model.rxnComps=cellfun(@toStr,raw(2:end,I(i)),'UniformOutput',false);
         case 10
-        	model.subSystems=cellfun(@toStr,raw(2:end,I(i)),'UniformOutput',false);
+        	subsystems=cellfun(@toStr,raw(2:end,I(i)),'UniformOutput',false);
         case 11
         	reactionReplacement=cellfun(@toStr,raw(2:end,I(i)),'UniformOutput',false);
         case 12
@@ -420,6 +420,9 @@
 if ~isempty(model.rxnConfidenceScores)
 	model.rxnConfidenceScores=str2double(model.rxnConfidenceScores);
 end
+for i=1:numel(subsystems)
+    model.subSystems{i,1}=cellstr(strsplit(subsystems{i,1},';'));
+end;
 
 %Check that all necessary reaction info has been loaded
 if isempty(equations)
diff --git a/readme/RAVEN_structure_fields.xlsx b/readme/RAVEN_structure_fields.xlsx
index a7c7e2f2..d51dc8b3 100644
Binary files a/readme/RAVEN_structure_fields.xlsx and b/readme/RAVEN_structure_fields.xlsx differ