Description of simplifyModel

0001 function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,...
0002     deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, constrainReversible, reservedRxns, suppressWarnings)
0003 % simplifyModel
0004 %   Simplifies a model by deleting reactions/metabolites
0005 %
0006 %   model                 a model structure
0007 %   deleteUnconstrained   delete metabolites marked as unconstrained (opt, default true)
0008 %   deleteDuplicates      delete all but one of duplicate reactions (opt, default false)
0009 %   deleteZeroInterval    delete reactions that are constrained to zero flux (opt, default false)
0010 %   deleteInaccessible    delete dead end reactions (opt, default false)
0011 %   deleteMinMax          delete reactions that cannot carry a flux by trying
0012 %                         to minimize/maximize the flux through that
0013 %                         reaction. May be time consuming (opt, default false)
0014 %   groupLinear           group linearly dependent pathways (opt, default false)
0015 %   constrainReversible   check if there are reversible reactions which can
0016 %                         only carry flux in one direction, and if so
0017 %                         constrain them to be irreversible. This tends to
0018 %                         allow for more reactions grouped when using
0019 %                         groupLinear (opt, default false)
0020 %   reservedRxns          cell array with reaction IDs that are not allowed to be
0021 %                         removed (opt)
0022 %   suppressWarnings      true if warnings should be suppressed (opt,
0023 %                         default false)
0024 %
0025 %   reducedModel          an updated model structure
0026 %   deletedReactions      a cell array with the IDs of all deleted reactions
0027 %   deletedMetabolites    a cell array with the IDs of all deleted
0028 %                         metabolites
0029 %
0030 %   This function is for reducing the model size by removing
0031 %   reactions and associated metabolites that cannot carry flux. It can also
0032 %   be used for identifying different types of gaps.
0033 %
0034 %   Usage: [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,...
0035 %           deleteUnconstrained, deleteDuplicates, deleteZeroInterval,...
0036 %           deleteInaccessible, deleteMinMax, groupLinear,...
0037 %           constrainReversible, reservedRxns, suppressWarnings)
0038 
0039 if nargin<2
0040     deleteUnconstrained=true;
0041 end
0042 if nargin<3
0043     deleteDuplicates=false;
0044 end
0045 if nargin<4
0046     deleteZeroInterval=false;
0047 end
0048 if nargin<5
0049     deleteInaccessible=false;
0050 end
0051 if nargin<6
0052     deleteMinMax=false;
0053 end
0054 if nargin<7
0055     groupLinear=false;
0056 end
0057 if nargin<8
0058     constrainReversible=false;
0059 end
0060 if nargin<9
0061     reservedRxns=[];
0062 else
0063     reservedRxns=convertCharArray(reservedRxns);
0064 end
0065 if nargin<10
0066     suppressWarnings=false;
0067 end
0068 
0069 reducedModel=model;
0070 deletedReactions={};
0071 deletedMetabolites={};
0072 
0073 if deleteUnconstrained==true
0074     if isfield(reducedModel,'unconstrained')
0075         %Remove unbalanced metabolites
0076         deletedMetabolites=reducedModel.mets(reducedModel.unconstrained~=0);
0077         reducedModel=removeMets(reducedModel,reducedModel.unconstrained~=0,false,false,false,true);
0078         reducedModel=rmfield(reducedModel,'unconstrained');
0079     end
0080 end
0081 
0082 if deleteDuplicates==true
0083     %Delete all but the last occurrence of duplicate reactions. The
0084     %reactions must have the same bounds, reversibility, and objective
0085     %coefficient to be regarded as duplicate
0086     [reducedModel, rxnsToDelete, ~]=contractModel(reducedModel);
0087     deletedReactions=[deletedReactions; rxnsToDelete];
0088 end
0089 
0090 if deleteZeroInterval==true
0091     %Find all reactions where both LB and UB are 0
0092     zeroIntervalReactions=and(reducedModel.lb==0,reducedModel.ub==0);
0093     
0094     rxnsToDelete=setdiff(reducedModel.rxns(zeroIntervalReactions),reservedRxns);
0095     deletedReactions=[deletedReactions; rxnsToDelete];
0096     
0097     %Remove reactions
0098     reducedModel=removeReactions(reducedModel,rxnsToDelete);
0099     
0100     %Find metabolites that no longer are used and delete them
0101     notInUse=sum(reducedModel.S~=0,2)==0;
0102     deletedMetabolites=[deletedMetabolites;reducedModel.mets(notInUse)];
0103     
0104     %Remove metabolites
0105     reducedModel=removeMets(reducedModel,notInUse);
0106 end
0107 
0108 if deleteInaccessible==true
0109     %Print a warning if exchange metabolites haven't been deleted yet. This
0110     %often means that the only allowed products are the metabolites that
0111     %are taken up be the system
0112     if isfield(reducedModel,'unconstrained') && suppressWarnings==false
0113         EM='Removing dead-end reactions before removing exchange metabolites';
0114         dispEM(EM,false);
0115     end
0116     
0117     while true
0118         %Find all metabolites that are only reactants or only products.
0119         %Delete those metabolites and reactions. This is done until no more
0120         %such reactions are found
0121         
0122         %Construct a stoichiometric matrix where all reactions are written
0123         %as reversible
0124         
0125         %Add fake exchange reactions if model.b~=0
0126         in=any(reducedModel.b<0,2);
0127         out=any(reducedModel.b>0,2);
0128         I=speye(numel(reducedModel.mets));
0129         revS=[reducedModel.S,reducedModel.S(:,reducedModel.lb<0)*-1 I(:,in) I(:,out)*-1];
0130         
0131         metUsage=sum(abs(revS')>0);
0132         onlyProducts=sum(revS'>0) == metUsage;
0133         onlyReactants=sum(revS'<0) == metUsage;
0134         
0135         %Also remove metabolites that only participate in one reversible
0136         %reaction Don't remove the ones in one reversible reaction if is
0137         %also has a non-zero coefficient in model.b
0138         notInUse=onlyProducts | onlyReactants | (sum(abs(reducedModel.S')>0)<=1 & (~in & ~out)');
0139         deletedRxn=false;
0140         if any(notInUse)
0141             %Find their corresponding reactions
0142             rxnsNotInUse=sum(abs(reducedModel.S(notInUse,:))>0,1)>0;
0143             rxnsToDelete=setdiff(reducedModel.rxns(rxnsNotInUse),reservedRxns);
0144             deletedReactions=[deletedReactions; rxnsToDelete];
0145             
0146             %Remove reactions
0147             reducedModel=removeReactions(reducedModel,rxnsToDelete);
0148             
0149             %Remove metabolites. Recalculate since it could be that some
0150             %cannot be deleted due to reserved rxns
0151             notInUse=sum(reducedModel.S~=0,2)==0;
0152             deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0153             reducedModel=removeMets(reducedModel,notInUse);
0154             
0155             %It could happen that it just finds reserved reactions and gets
0156             %stuck
0157             if ~isempty(rxnsToDelete)
0158                 deletedRxn=true;
0159             end
0160         else
0161             break;
0162         end
0163         if deletedRxn==false
0164             break;
0165         end
0166     end
0167 end
0168 
0169 if deleteMinMax==true
0170     %Get reactions that can't carry fluxes. This should be done
0171     %algebraically if possible
0172     I=~haveFlux(reducedModel);
0173     
0174     %Remove reactions
0175     rxnsToDelete=setdiff(reducedModel.rxns(I),reservedRxns);
0176     deletedReactions=[deletedReactions; rxnsToDelete];
0177     reducedModel=removeReactions(reducedModel,rxnsToDelete);
0178     
0179     %Remove metabolites
0180     notInUse=sum(reducedModel.S~=0,2)==0;
0181     deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0182     reducedModel=removeMets(reducedModel,notInUse);
0183 end
0184 
0185 if constrainReversible==true
0186     revs=find(reducedModel.rev);
0187     [I,J]=getAllowedBounds(reducedModel,revs);
0188     
0189     I=abs(I);
0190     J=abs(J);
0191     
0192     %Get the "small" values
0193     K=I<10^-10;
0194     L=J<10^-10;
0195     
0196     %Keep the values where only one direction is zero (small)
0197     I(K==L)=[];
0198     J(K==L)=[];
0199     revs(K==L)=[];
0200     
0201     %Change the reversibility of the remaining reactions
0202     reducedModel.rev(revs(J>10^-10))=0; %Ignore reverse direction
0203     reducedModel.lb(revs(J>10^-10))=0;
0204     
0205     toSwitch=revs(I>10^-10);
0206     reducedModel.rev(toSwitch)=0; %Change directionality
0207     reducedModel.ub(toSwitch)=reducedModel.lb(toSwitch)*-1;
0208     reducedModel.lb(toSwitch)=0;
0209     reducedModel.S(:,toSwitch)=reducedModel.S(:,toSwitch).*-1;
0210     reducedModel.c(toSwitch)=reducedModel.c(toSwitch)*-1;
0211 end
0212 if groupLinear==true
0213     if  suppressWarnings==false
0214         fprintf('NOTE: You have chosen to group linear reactions. This option does not keep track of gene/reaction associations when reactions are merged. Deleting all gene information\n');
0215     end
0216     
0217     reducedModel.genes={};
0218     reducedModel.rxnGeneMat=sparse(numel(reducedModel.rxns),0);
0219     reducedModel.grRules(:)={''};
0220     
0221     if isfield(reducedModel,'geneShortNames')
0222         reducedModel.geneShortNames={};
0223     end
0224     if isfield(reducedModel,'geneMiriams')
0225         reducedModel.geneMiriams={};
0226     end
0227     if isfield(reducedModel,'geneComps')
0228         reducedModel.geneComps=[];
0229     end
0230     
0231     %Convert the model to irreversible
0232     irrevModel=convertToIrrev(reducedModel);
0233     
0234     %Loop through and iteratively group linear reactions
0235     while 1
0236         %Get the banned reaction indexes. Note that the indexes will change
0237         %in each iteration, but the names will not as they won't be merged
0238         %with any other reaction
0239         bannedIndexes=getIndexes(irrevModel,reservedRxns,'rxns');
0240         
0241         %Select all metabolites that are only present as reactants/products
0242         %in one reaction
0243         singleNegative=find(sum(irrevModel.S'<0)==1);
0244         singlePositive=find(sum(irrevModel.S'>0)==1);
0245         
0246         %Retrieve the common metabolites
0247         common=intersect(singleNegative,singlePositive);
0248         
0249         mergedSome=false;
0250         
0251         %Loop through each of them and see if the reactions should be
0252         %merged
0253         for i=1:numel(common)
0254             involvedRxns=find(irrevModel.S(common(i),:));
0255             
0256             %Check so that one or both of the reactions haven't been merged
0257             %already
0258             if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns))
0259                 %Calculate how many times the second reaction has to be
0260                 %multiplied before being merged with the first
0261                 stoichRatio=abs(irrevModel.S(common(i),involvedRxns(1))/irrevModel.S(common(i),involvedRxns(2)));
0262                 
0263                 %Add the second to the first
0264                 irrevModel.S(:,involvedRxns(1))=irrevModel.S(:,involvedRxns(1))+irrevModel.S(:,involvedRxns(2))*stoichRatio;
0265                 
0266                 %Clear the second reaction
0267                 irrevModel.S(:,involvedRxns(2))=0;
0268                 
0269                 %This is to prevent numerical issues. It should be 0
0270                 %already
0271                 irrevModel.S(common(i),involvedRxns(1))=0;
0272                 
0273                 %At this point the second reaction is certain to be deleted
0274                 %in a later step and can therefore be ignored
0275                 
0276                 %Recalculate the bounds for the new reaction. This can be
0277                 %problematic since the scale of the bounds may change
0278                 %dramatically. Let the most constraining reaction determine
0279                 %the new bound
0280                 lb1=irrevModel.lb(involvedRxns(1));
0281                 lb2=irrevModel.lb(involvedRxns(2));
0282                 ub1=irrevModel.ub(involvedRxns(1));
0283                 ub2=irrevModel.ub(involvedRxns(2));
0284                 
0285                 if lb2~=-inf
0286                     irrevModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio);
0287                 end
0288                 if ub2~=inf
0289                     irrevModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio);
0290                 end
0291                 
0292                 %Then recalculate the objective coefficient. The resulting
0293                 %coefficient is the weighted sum of the previous
0294                 irrevModel.c(involvedRxns(1))=irrevModel.c(involvedRxns(1))+irrevModel.c(involvedRxns(2))*stoichRatio;
0295                 
0296                 %Iterate again
0297                 mergedSome=true;
0298             end
0299         end
0300         
0301         %All possible reactions merged
0302         if mergedSome==false
0303             break;
0304         end
0305         
0306         %Now delete all reactions that involve no metabolites
0307         I=find(sum(irrevModel.S~=0)==0);
0308         
0309         %Remove reactions
0310         irrevModel=removeReactions(irrevModel,I);
0311         
0312         %Remove metabolites
0313         notInUse=sum(irrevModel.S~=0,2)==0;
0314         irrevModel=removeMets(irrevModel,notInUse);
0315     end
0316     
0317     reducedModel=irrevModel;
0318 end
0319 end
simplifyModel

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
