Home > RAVEN > getElementalBalance.m

getElementalBalance

PURPOSE ^

getElementalBalance

SYNOPSIS ^

function balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)

DESCRIPTION ^

 getElementalBalance
   Checks a model to see if the reactions are elementally balanced

   model             a model structure
   rxns              either a cell array of reaction IDs, a logical vector 
                     with the same number of elements as reactions in the model,
                     of a vector of indexes. Only these reactions will be
                     checked (opt, default model.rxns)
   printUnbalanced   print warnings about the reactions that were 
                     unbalanced (opt, default false)
   printUnparsable   print warnings about the reactions that cannot be
                     parsed (opt, default false)

   balanceStructure
       balanceStatus    1 if the reaction is balanced, 0 if it's unbalanced, 
                      -1 if it couldn't be balanced due to missing information, 
                      -2 if it couldn't be balanced due to an error
       elements
           abbrevs     cell array with abbreviations for all used elements
           names       cell array with the names for all used elements
       leftComp        MxN matrix with the sum of coefficients for each of
                       the elements (N) for the left side of the
                       reactions (M)
       rightComp       the corresponding matrix for the right side

   Usage: balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)

   Rasmus Agren, 2013-08-01

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)
0002 % getElementalBalance
0003 %   Checks a model to see if the reactions are elementally balanced
0004 %
0005 %   model             a model structure
0006 %   rxns              either a cell array of reaction IDs, a logical vector
0007 %                     with the same number of elements as reactions in the model,
0008 %                     of a vector of indexes. Only these reactions will be
0009 %                     checked (opt, default model.rxns)
0010 %   printUnbalanced   print warnings about the reactions that were
0011 %                     unbalanced (opt, default false)
0012 %   printUnparsable   print warnings about the reactions that cannot be
0013 %                     parsed (opt, default false)
0014 %
0015 %   balanceStructure
0016 %       balanceStatus    1 if the reaction is balanced, 0 if it's unbalanced,
0017 %                      -1 if it couldn't be balanced due to missing information,
0018 %                      -2 if it couldn't be balanced due to an error
0019 %       elements
0020 %           abbrevs     cell array with abbreviations for all used elements
0021 %           names       cell array with the names for all used elements
0022 %       leftComp        MxN matrix with the sum of coefficients for each of
0023 %                       the elements (N) for the left side of the
0024 %                       reactions (M)
0025 %       rightComp       the corresponding matrix for the right side
0026 %
0027 %   Usage: balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)
0028 %
0029 %   Rasmus Agren, 2013-08-01
0030 %
0031 
0032 if nargin<2
0033     rxns=[];
0034 end
0035 
0036 if nargin<3
0037     printUnbalanced=false;
0038 end
0039 
0040 if nargin<4
0041     printUnparsable=false;
0042 end
0043 
0044 if ~isempty(rxns)
0045    indexes=~getIndexes(model,rxns,'rxns',true);
0046    model=removeRxns(model,indexes,true);
0047 end
0048 
0049 balanceStructure.balanceStatus=nan(numel(model.rxns),1);
0050 
0051 %Get the formulas
0052 if isfield(model,'metFormulas')
0053     [balanceStructure.elements, useMat, exitFlag]=parseFormulas(model.metFormulas, true);
0054 else
0055     if isfield(model,'inchis')
0056         [balanceStructure.elements, useMat, exitFlag]=parseFormulas(model.inchis, true,true);
0057     else
0058         dispEM('The model must contain either the "metFormulas" or the "inchis" field in order to test for elemental balancing');
0059     end
0060 end
0061 
0062 balanceStructure.leftComp=zeros(numel(model.rxns),numel(balanceStructure.elements.names));
0063 balanceStructure.rightComp=balanceStructure.leftComp;
0064 
0065 %Look at the right and left side of the reactions separately
0066 S{1}=model.S;
0067 S{2}=model.S;
0068 S{1}(S{1}>0)=0; %Left side
0069 S{2}(S{2}<0)=0; %Right side
0070 S{1}=abs(S{1}); %Both should have positive values
0071 
0072 %First do it for left side and then for right
0073 for i=1:2
0074     for j=1:numel(model.rxns)
0075         %Get the balance status of the involved mets
0076         I=exitFlag(S{i}(:,j)~=0);
0077         if any(I==-1)
0078             balanceStructure.balanceStatus(j)=-2;
0079         end
0080         if any(I==0)
0081             %Don't change a more serious error to a less serious one
0082             balanceStructure.balanceStatus(j)=min(-1,balanceStructure.balanceStatus(j));
0083         end
0084         %Loop through each element
0085         for k=1:numel(balanceStructure.elements.names)
0086             if i==1
0087                 balanceStructure.leftComp(j,k)=sum(S{i}(:,j).*useMat(:,k));
0088             else
0089                 balanceStructure.rightComp(j,k)=sum(S{i}(:,j).*useMat(:,k));
0090             end
0091         end
0092     end
0093 end
0094 
0095 %Now compare the left and right sides to find which are unbalanced. This is
0096 %done even if the reaction as a whole couldn't be balanced
0097 total=abs(balanceStructure.rightComp-balanceStructure.leftComp)>10^-8; %To deal with roundoff error
0098 
0099 %Get which reactions are unbalanced. Don't change an error to just
0100 %unbalanced
0101 balanceStructure.balanceStatus(any(total,2))=min(balanceStructure.balanceStatus(any(total,2)),0);
0102 
0103 %The remaining ones are all balanced
0104 balanceStructure.balanceStatus(isnan(balanceStructure.balanceStatus))=1;
0105 
0106 %Print warnings
0107 toPrint=[];
0108 if printUnbalanced==true
0109    toPrint=[toPrint;find(balanceStructure.balanceStatus==0)];
0110 end
0111 if printUnparsable==true
0112    toPrint=[toPrint;find(balanceStructure.balanceStatus<0)];
0113 end
0114 
0115 toPrint=sort(toPrint);
0116 for i=1:numel(toPrint)
0117     if balanceStructure.balanceStatus(toPrint(i))<0
0118         if balanceStructure.balanceStatus(toPrint(i))==-1
0119             dispEM(['The reaction ' model.rxns{toPrint(i)} ' could not be balanced due to missing information'],false);
0120         else
0121             dispEM(['The reaction ' model.rxns{toPrint(i)} ' could not be balanced due to a parsing error'],false);
0122         end
0123     else
0124        %Find the compounds that it's not balanced for
0125        notBalanced=find(total(toPrint(i),:));
0126        for j=1:numel(notBalanced)
0127             dispEM(['The reaction ' model.rxns{toPrint(i)} ' is not balanced with respect to ' balanceStructure.elements.names{notBalanced(j)}],false);
0128        end
0129     end
0130 end
0131 end

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005