0001 function [massImbalance,imBalancedMass,imBalancedCharge,imBalancedBool,Elements] = checkMassChargeBalance(model,rxnBool,printLevel)
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
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038 [nMet,nRxn]=size(model.S);
0039 if exist('rxnBool','var')
0040 if ~isempty(rxnBool)
0041 if length(rxnBool)~=nRxn
0042 rxnBool2=false(nRxn,1);
0043 rxnBool2(rxnBool)=1;
0044 rxnBool=rxnBool2;
0045 end
0046 else
0047 model=findSExRxnInd(model);
0048
0049 rxnBool=model.SIntRxnBool;
0050 end
0051 else
0052 model=findSExRxnInd(model);
0053
0054 rxnBool=model.SIntRxnBool;
0055 end
0056 if ~exist('printLevel','var')
0057 printLevel=0;
0058 end
0059
0060
0061 Elements = {'H','C', 'O', 'P', 'S', 'N', 'Mg','X','Fe','Zn','Co','R'};
0062
0063 E=sparse(nMet,length(Elements));
0064 massImbalance=sparse(nRxn,length(Elements));
0065 for j = 1 : length(Elements)
0066 if j==1
0067 [dE,E_el]=checkBalance(model,Elements{j},printLevel);
0068 massImbalance(:,j)=dE;
0069 E(:,j)=E_el;
0070 fprintf('%s\n',['Checked element ' Elements{j}]);
0071 else
0072
0073
0074 [massImbalance(:,j),E(:,j)]=checkBalance(model,Elements{j},0);
0075 fprintf('%s\n',['Checking element ' Elements{j}]);
0076 end
0077 end
0078 massImbalance(~rxnBool,:)=0;
0079 imBalancedBool=sum(abs(massImbalance'))'~=0;
0080
0081 imBalancedBool=rxnBool & imBalancedBool;
0082
0083 imBalancedMass=cell(nRxn,1);
0084 for i = 1 : nRxn
0085 imBalancedMass{i,1}='';
0086 if imBalancedBool(i)
0087 for j = 1 : length(Elements)
0088 if massImbalance(i,j)~=0
0089 if ~strcmp(imBalancedMass{i,1},'')
0090 imBalancedMass{i,1} = [imBalancedMass{i,1} ', ' int2str(massImbalance(i,j)) ' ' Elements{j}];
0091 else
0092 imBalancedMass{i,1} = [int2str(massImbalance(i,j)) ' ' Elements{j}];
0093 end
0094 end
0095
0096 end
0097 if strfind(imBalancedMass{i,1},'NaN')
0098 imBalancedMass{i,1}='NaN';
0099 end
0100 end
0101 if mod(i,1000)==0
0102 fprintf('%n\t%s\n',i,['reactions checked for ' Elements{j} ' balance']);
0103 end
0104 end
0105 if printLevel==-1
0106 firstMissing=0;
0107 for p=1:nRxn
0108 if ~strcmp(imBalancedMass{p,1},'')
0109
0110
0111 if ~strcmp(imBalancedMass{p,1},'NaN')
0112 if ~firstMissing
0113 fid=fopen('mass_imbalanced_reactions.txt','w');
0114 fprintf(fid,'%s;%s;%s;%s\n','#Rxn','rxnAbbr','imbalance','equation');
0115
0116 warning('There are mass imbalanced reactions, see mass_imbalanced_reactions.txt')
0117 firstMissing=1;
0118 end
0119 equation=printRxnFormula(model,model.rxns(p),0);
0120 fprintf(fid,'%s;%s;%s;%s\n',int2str(p),model.rxns{p},imBalancedMass{p,1},equation{1});
0121 for m=1:size(model.S,1)
0122 if model.S(m,p)~=0
0123 fprintf(fid,'%s\t%s\t%s\t%s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,p)),int2str(E(m)),model.metFormulas{m});
0124 end
0125 end
0126 end
0127 end
0128 end
0129 if firstMissing
0130 fclose(fid);
0131 end
0132 end
0133 if printLevel==1
0134 for p=1:nRxn
0135 if ~strcmp(imBalancedMass{p,1},'')
0136
0137
0138 if ~strcmp(imBalancedMass{p,1},'NaN')
0139 equation=printRxnFormula(model,model.rxns(p),0);
0140 fprintf('%6s\t%30s\t%10s\t%s\n',int2str(p),model.rxns{p},imBalancedMass{p,1},equation{1});
0141 if 0
0142 for m=1:size(model.S,1)
0143 if model.S(m,p)~=0
0144 fprintf(fid,'%s\t%s\t%s\t%s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,p)),int2str(E(m)),model.metFormulas{m});
0145 end
0146 end
0147 end
0148 end
0149 end
0150 end
0151 end
0152
0153
0154 if nnz(strcmp('',imBalancedMass))==nRxn
0155 imBalancedMass=[];
0156 end
0157
0158
0159 imBalancedCharge=[];
0160 firstMissing=0;
0161 if isfield(model, 'metCharges')
0162 for m=1:nMet
0163 if isnan(model.metCharges(m)) && ~isempty(model.metFormulas{m})
0164 if printLevel==1
0165 fprintf('%s\t%s\n',int2str(m),[model.mets{m} ' has no charge but has formula.'])
0166 if ~firstMissing
0167 warning('model structure must contain model.metCharges field for each metabolite');
0168 end
0169 firstMissing=1;
0170 end
0171 if printLevel==-1
0172 if ~firstMissing
0173 fid=fopen('metabolites_without_charge.txt','w');
0174 end
0175 firstMissing=1;
0176 fprintf(fid,'%s\t%s\n',int2str(m),model.mets{m})
0177 end
0178 else
0179 dC=model.S'*model.metCharges;
0180 end
0181 end
0182 if any(dC(rxnBool))~=0
0183 imBalancedCharge=dC;
0184 imBalancedCharge(~rxnBool)=0;
0185 else
0186 imBalancedCharge=[];
0187 end
0188 end
0189
0190 if printLevel==-1
0191 firstMissing=0;
0192 if ~isempty(imBalancedCharge)
0193 for q=1:nRxn
0194 if model.SIntRxnBool(q) && dC(q)~=0 && strcmp(imBalancedMass{p,1},'')
0195 if ~firstMissing
0196 fid=fopen('charge_imbalanced_reactions.txt','w');
0197 warning('There are charged imbalanced reactions (that are mass balanced), see charge_imbalanced_reactions.txt')
0198 firstMissing=1;
0199 end
0200 equation=printRxnFormula(model,model.rxns(q),0);
0201 fprintf(fid,'%s\t%s\t%s\n',int2str(q),model.rxns{q},equation{1});
0202 if 0
0203 for m=1:size(model.S,1)
0204 if model.S(m,q)~=0
0205 fprintf(fid,'%s\t%15s\t%3s\t%3s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,q)),int2str(model.metCharges(m)),model.metFormulas{m});
0206 end
0207 end
0208 end
0209 end
0210 end
0211 if firstMissing
0212 fclose(fid);
0213 end
0214 end
0215 end
0216
0217 if printLevel==1
0218 if ~isempty(imBalancedCharge)
0219 fprintf('%s\n','Mass balanced, but charged imbalanced reactions:')
0220 for q=1:nRxn
0221 if model.SIntRxnBool(q) && dC(q)~=0 && strcmp(imBalancedMass{p,1},'')
0222 equation=printRxnFormula(model,model.rxns(q),0);
0223 fprintf('%s\t%s\t%s\n',int2str(q),model.rxns{q},equation{1});
0224 if 1
0225 for m=1:size(model.S,1)
0226 if model.S(m,q)~=0
0227 fprintf('%s\t%15s\t%3s\t%3s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,q)),int2str(model.metCharges(m)),model.metFormulas{m});
0228 end
0229 end
0230 end
0231 end
0232 end
0233 end
0234 end
0235
0236 if ~isempty(imBalancedCharge)
0237 imBalancedBool = imBalancedBool | imBalancedCharge~=0;
0238 end
0239
0240
0241
0242
0243