getElementalBalance

 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 (optional, default model.rxns)
   printUnbalanced   print warnings about the reactions that were
                     unbalanced (optional, default false)
   printUnparsable   print warnings about the reactions that cannot be
                     parsed (optional, 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)

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 (optional, default model.rxns)
0010 %   printUnbalanced   print warnings about the reactions that were
0011 %                     unbalanced (optional, default false)
0012 %   printUnparsable   print warnings about the reactions that cannot be
0013 %                     parsed (optional, 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 if nargin<2
0030     rxns=[];
0031 elseif ~islogical(rxns) && ~isnumeric(rxns)
0032     rxns=convertCharArray(rxns);
0033 end
0034 
0035 if nargin<3
0036     printUnbalanced=false;
0037 end
0038 
0039 if nargin<4
0040     printUnparsable=false;
0041 end
0042 
0043 if ~isempty(rxns)
0044     indexes=~getIndexes(model,rxns,'rxns',true);
0045     model=removeReactions(model,indexes,true);
0046 end
0047 
0048 balanceStructure.balanceStatus=nan(numel(model.rxns),1);
0049 
0050 %Get the formulas
0051 if isfield(model,'metFormulas')
0052     [balanceStructure.elements, useMat, exitFlag]=parseFormulas(model.metFormulas, true);
0053 else
0054     if isfield(model,'inchis')
0055         [balanceStructure.elements, useMat, exitFlag]=parseFormulas(model.inchis, true,true);
0056     else
0057         EM='The model must contain either the "metFormulas" or the "inchis" field in order to test for elemental balancing';
0058         dispEM(EM);
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             EM=['The reaction ' model.rxns{toPrint(i)} ' could not be balanced due to missing information'];
0120             dispEM(EM,false);
0121         else
0122             EM=['The reaction ' model.rxns{toPrint(i)} ' could not be balanced due to a parsing error'];
0123             dispEM(EM,false);
0124         end
0125     else
0126         %Find the compounds that it's not balanced for
0127         notBalanced=find(total(toPrint(i),:));
0128         for j=1:numel(notBalanced)
0129             EM=['The reaction ' model.rxns{toPrint(i)} ' is not balanced with respect to ' balanceStructure.elements.names{notBalanced(j)}];
0130             dispEM(EM,false);
0131         end
0132     end
0133 end
0134 
0135 % Re-order the structure entries so they're consistent with the ordering of
0136 % the input reaction indexes
0137 if ~isempty(rxns)
0138     rxns = getIndexes(model,rxns,'rxns');
0139     [~,i] = sort(rxns);
0140     balanceStructure.balanceStatus(i) = balanceStructure.balanceStatus;
0141     balanceStructure.leftComp(i,:) = balanceStructure.leftComp;
0142     balanceStructure.rightComp(i,:) = balanceStructure.rightComp;
0143 end