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 (optional, default true)
0008 %   deleteDuplicates      delete all but one of duplicate reactions (optional, default false)
0009 %   deleteZeroInterval    delete reactions that are constrained to zero flux (optional, default false)
0010 %   deleteInaccessible    delete dead end reactions (optional, 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 (optional, default false)
0014 %   groupLinear           group linearly dependent pathways (optional, 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 (optional, default false)
0020 %   reservedRxns          cell array with reaction IDs that are not allowed to be
0021 %                         removed (optional)
0022 %   suppressWarnings      true if warnings should be suppressed (optional,
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,'proteins')
0225         reducedModel.proteins={};
0226     end
0227     if isfield(reducedModel,'geneMiriams')
0228         reducedModel.geneMiriams={};
0229     end
0230     if isfield(reducedModel,'geneComps')
0231         reducedModel.geneComps=[];
0232     end
0233     
0234     %Convert the model to irreversible
0235     irrevModel=convertToIrrev(reducedModel);
0236     
0237     %Loop through and iteratively group linear reactions
0238     while 1
0239         %Get the banned reaction indexes. Note that the indexes will change
0240         %in each iteration, but the names will not as they won't be merged
0241         %with any other reaction
0242         bannedIndexes=getIndexes(irrevModel,reservedRxns,'rxns');
0243         
0244         %Select all metabolites that are only present as reactants/products
0245         %in one reaction
0246         singleNegative=find(sum(irrevModel.S'<0)==1);
0247         singlePositive=find(sum(irrevModel.S'>0)==1);
0248         
0249         %Retrieve the common metabolites
0250         common=intersect(singleNegative,singlePositive);
0251         
0252         mergedSome=false;
0253         
0254         %Loop through each of them and see if the reactions should be
0255         %merged
0256         for i=1:numel(common)
0257             involvedRxns=find(irrevModel.S(common(i),:));
0258             
0259             %Check so that one or both of the reactions haven't been merged
0260             %already
0261             if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns))
0262                 %Calculate how many times the second reaction has to be
0263                 %multiplied before being merged with the first
0264                 stoichRatio=abs(irrevModel.S(common(i),involvedRxns(1))/irrevModel.S(common(i),involvedRxns(2)));
0265                 
0266                 %Add the second to the first
0267                 irrevModel.S(:,involvedRxns(1))=irrevModel.S(:,involvedRxns(1))+irrevModel.S(:,involvedRxns(2))*stoichRatio;
0268                 
0269                 %Clear the second reaction
0270                 irrevModel.S(:,involvedRxns(2))=0;
0271                 
0272                 %This is to prevent numerical issues. It should be 0
0273                 %already
0274                 irrevModel.S(common(i),involvedRxns(1))=0;
0275                 
0276                 %At this point the second reaction is certain to be deleted
0277                 %in a later step and can therefore be ignored
0278                 
0279                 %Recalculate the bounds for the new reaction. This can be
0280                 %problematic since the scale of the bounds may change
0281                 %dramatically. Let the most constraining reaction determine
0282                 %the new bound
0283                 lb1=irrevModel.lb(involvedRxns(1));
0284                 lb2=irrevModel.lb(involvedRxns(2));
0285                 ub1=irrevModel.ub(involvedRxns(1));
0286                 ub2=irrevModel.ub(involvedRxns(2));
0287                 
0288                 if lb2~=-inf
0289                     irrevModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio);
0290                 end
0291                 if ub2~=inf
0292                     irrevModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio);
0293                 end
0294                 
0295                 %Then recalculate the objective coefficient. The resulting
0296                 %coefficient is the weighted sum of the previous
0297                 irrevModel.c(involvedRxns(1))=irrevModel.c(involvedRxns(1))+irrevModel.c(involvedRxns(2))*stoichRatio;
0298                 
0299                 %Iterate again
0300                 mergedSome=true;
0301             end
0302         end
0303         
0304         %All possible reactions merged
0305         if mergedSome==false
0306             break;
0307         end
0308         
0309         %Now delete all reactions that involve no metabolites
0310         I=find(sum(irrevModel.S~=0)==0);
0311         
0312         %Remove reactions
0313         irrevModel=removeReactions(irrevModel,I);
0314         
0315         %Remove metabolites
0316         notInUse=sum(irrevModel.S~=0,2)==0;
0317         irrevModel=removeMets(irrevModel,notInUse);
0318     end
0319     
0320     reducedModel=irrevModel;
0321 end
0322 end
simplifyModel

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
