getAllSubGraphs

 getAllSubGraphs
   Get all metabolic subgraphs in a model. Two metabolites 
   are connected if they share a reaction.

   Input:
   model         a model structure

   Output:
   subGraphs     a boolean matrix where the rows correspond to the metabolites
                 and the columns to which subgraph they are assigned to. The
                 columns are ordered so that larger subgraphs come first

 Usage: subGraphs=getAllSubGraphs(model)

0001 function subGraphs=getAllSubGraphs(model)
0002 % getAllSubGraphs
0003 %   Get all metabolic subgraphs in a model. Two metabolites
0004 %   are connected if they share a reaction.
0005 %
0006 %   Input:
0007 %   model         a model structure
0008 %
0009 %   Output:
0010 %   subGraphs     a boolean matrix where the rows correspond to the metabolites
0011 %                 and the columns to which subgraph they are assigned to. The
0012 %                 columns are ordered so that larger subgraphs come first
0013 %
0014 % Usage: subGraphs=getAllSubGraphs(model)
0015 
0016 %Generate the connectivity graph. Metabolites are connected through
0017 %reactions. This is not a bipartite graph with the reactions.
0018 G=model.S*model.S';
0019 G(G~=0)=1;
0020 
0021 %Keeps track of which mets have been assigned to a subgraph
0022 isAssigned=false(numel(model.mets),1);
0023 
0024 %Allocate space for subgraphs, initially one graph for each met
0025 subGraphs=false(numel(model.mets),numel(model.mets));
0026 
0027 %Main loop continues until all mets have been assigned to a subgraph
0028 counter=1;
0029 while ~all(isAssigned)
0030     currentSG=false(numel(model.mets),1);
0031     %Find the first non-assigned metabolite and assign it to the current SG
0032     I=find(isAssigned==false,1);
0033     currentSG(I)=true;
0034     isAssigned(I)=true;
0035     %Then iteratively add all mets that are connected to each other, until
0036     %no more such mets can be found
0037     while true
0038         J=sum(G(:,currentSG),2);
0039         
0040         %Get the new mets for this SG. Also assign them to the current SG
0041         newOnes=J~=0 & currentSG==false;
0042         isAssigned(newOnes)=true;
0043         currentSG(newOnes)=true;
0044         
0045         %If there are no new mets, then abort
0046         if ~any(newOnes)
0047             subGraphs(currentSG,counter)=true;
0048             counter=counter+1;
0049             break;
0050         end
0051     end
0052 end
0053 subGraphs=subGraphs(:,1:counter);
0054 
0055 [~, I]=sort(sum(subGraphs),'descend');
0056 
0057 subGraphs=subGraphs(:,I);
0058 
0059 %Also remove empty subgraphs (can happen when metabolites are never used)
0060 subGraphs(:,sum(subGraphs)==0)=[];
0061 end