0001 function balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 if nargin<2
0029 rxns=[];
0030 end
0031
0032 if nargin<3
0033 printUnbalanced=false;
0034 end
0035
0036 if nargin<4
0037 printUnparsable=false;
0038 end
0039
0040 if ~isempty(rxns)
0041 indexes=~getIndexes(model,rxns,'rxns',true);
0042 model=removeRxns(model,indexes);
0043 end
0044
0045
0046
0047 if isfield(model,'metFormulas')
0048
0049 [nC, nP, nS, nN, nR, nO, nH, exitFlag]=parseFormulas(model.metFormulas,true);
0050 else
0051 nC=zeros(numel(model.mets),1);
0052 nP=nC;
0053 nS=nC;
0054 nN=nC;
0055 nR=nC;
0056 nO=nC;
0057 nH=nC;
0058 exitFlag=nC;
0059 end
0060 if isfield(model,'inchis')
0061 [nCI, nPI, nSI, nNI, nRI, nOI, nHI, exitFlagI]=parseFormulas(model.inchis,true,true);
0062
0063
0064 nC(exitFlagI==1)=nCI(exitFlagI==1);
0065 nP(exitFlagI==1)=nPI(exitFlagI==1);
0066 nS(exitFlagI==1)=nSI(exitFlagI==1);
0067 nN(exitFlagI==1)=nNI(exitFlagI==1);
0068 nR(exitFlagI==1)=nRI(exitFlagI==1);
0069 nO(exitFlagI==1)=nOI(exitFlagI==1);
0070 nH(exitFlagI==1)=nHI(exitFlagI==1);
0071 exitFlag(exitFlagI==1)=exitFlagI(exitFlagI==1);
0072 end
0073
0074 balanceStructure.balanceStatus=zeros(numel(model.rxns),1);
0075 balanceStructure.leftComp=zeros(numel(model.rxns),7);
0076 balanceStructure.rightComp=zeros(numel(model.rxns),7);
0077
0078
0079 for i=1:numel(model.rxns)
0080 leftMets=find(model.S(:,i)<0);
0081 rightMets=find(model.S(:,i)>0);
0082
0083
0084
0085 balanceStructure.leftComp(i,1)=abs(sum(model.S(leftMets,i).*nC(leftMets)));
0086 balanceStructure.leftComp(i,2)=abs(sum(model.S(leftMets,i).*nP(leftMets)));
0087 balanceStructure.leftComp(i,3)=abs(sum(model.S(leftMets,i).*nS(leftMets)));
0088 balanceStructure.leftComp(i,4)=abs(sum(model.S(leftMets,i).*nN(leftMets)));
0089 balanceStructure.leftComp(i,5)=abs(sum(model.S(leftMets,i).*nR(leftMets)));
0090 balanceStructure.leftComp(i,6)=abs(sum(model.S(leftMets,i).*nO(leftMets)));
0091 balanceStructure.leftComp(i,7)=abs(sum(model.S(leftMets,i).*nH(leftMets)));
0092
0093 balanceStructure.rightComp(i,1)=abs(sum(model.S(rightMets,i).*nC(rightMets)));
0094 balanceStructure.rightComp(i,2)=abs(sum(model.S(rightMets,i).*nP(rightMets)));
0095 balanceStructure.rightComp(i,3)=abs(sum(model.S(rightMets,i).*nS(rightMets)));
0096 balanceStructure.rightComp(i,4)=abs(sum(model.S(rightMets,i).*nN(rightMets)));
0097 balanceStructure.rightComp(i,5)=abs(sum(model.S(rightMets,i).*nR(rightMets)));
0098 balanceStructure.rightComp(i,6)=abs(sum(model.S(rightMets,i).*nO(rightMets)));
0099 balanceStructure.rightComp(i,7)=abs(sum(model.S(rightMets,i).*nH(rightMets)));
0100
0101
0102 if all(exitFlag([leftMets; rightMets]))
0103 if all(exitFlag([leftMets; rightMets])==1)
0104
0105
0106 if sum(abs(balanceStructure.leftComp(i,:)-balanceStructure.rightComp(i,:)))<10^-6
0107 balanceStructure.balanceStatus(i)=1;
0108 end
0109 else
0110 balanceStructure.balanceStatus(i)=-2;
0111 end
0112 else
0113 balanceStructure.balanceStatus(i)=-1;
0114 end
0115 end
0116
0117
0118 toPrint=[];
0119 if printUnbalanced==true
0120 toPrint=[toPrint;find(balanceStructure.balanceStatus==0)];
0121 end
0122 if printUnparsable==true
0123 toPrint=[toPrint;find(balanceStructure.balanceStatus<0)];
0124 end
0125
0126 compounds={'carbon' 'phosphate' 'sulfate' 'nitrogen' 'R-groups' 'oxygen' 'hydrogen'};
0127 toPrint=sort(toPrint);
0128 for i=1:numel(toPrint)
0129 if balanceStructure.balanceStatus(toPrint(i))<0
0130 fprintf(['WARNING: The reaction ' model.rxns{toPrint(i)} ' could not be parsed\n']);
0131 else
0132
0133 notBalanced=find(abs(balanceStructure.leftComp(toPrint(i),:)-balanceStructure.rightComp(toPrint(i),:))>10^-6);
0134 for j=1:numel(notBalanced)
0135 fprintf(['WARNING: The reaction ' model.rxns{toPrint(i)} ' is not balanced with respect to ' compounds{notBalanced(j)} '\n']);
0136 end
0137 end
0138 end
0139
0140 end