Mass balance protons for each reaction by adjusting the hydrogen ion stoichiometric coefficient. For each non transport reaction the proton stoichiometric coefficient for each reaction is given by the difference between the sum of the average number of protons bound by substrate reactants and the average number of protons bound by product reactants i.e. proton S_ij = sum(substrates(S_ij*aveHbound)) - sum(products(S_ij*aveHbound)) For each transport reaction, model transport across the membrane as three reactions, one where non-proton reactants convert into non-proton metabolites involved in one compartment,then in the other compartment, non-proton metabolites convert into non-proton reactants. In between is a reconstruction transport reaction which must be elementally balanced for H to begin with. We assume that the transporter involved has affinity only for a particular species of metabolite defined by the reconstrution. INPUT model.S model.mets model.SIntRxnBool Boolean of internal reactions model.met(m).aveHbound average number of bound hydrogen ions model.met(m).formula average number of bound hydrogen ions model.met(m).charge charge of metabolite species from reconstruction OPTIONAL INPUT massImbalance nRxn x nElement matrix where massImbalance(i,j) is imbalance of reaction i for element j. massImbalance(i,j)==0 if balanced for that element. The first column is assumed to correspond to H balancing. OUTPUT model.S Stoichiometric matrix balanced for protons where each row. corresponds to a reactant at specific pH. Ronan M. T. Fleming
0001 function model=pHbalanceProtons(model,massImbalance) 0002 % Mass balance protons for each reaction by adjusting the hydrogen ion stoichiometric coefficient. 0003 % 0004 % For each non transport reaction the proton stoichiometric coefficient for each 0005 % reaction is given by the difference between the sum of the average number 0006 % of protons bound by substrate reactants and the average number of protons 0007 % bound by product reactants 0008 % 0009 % i.e. proton S_ij = sum(substrates(S_ij*aveHbound)) - 0010 % sum(products(S_ij*aveHbound)) 0011 % 0012 % For each transport reaction, model transport across the membrane as three 0013 % reactions, one where non-proton reactants convert into non-proton metabolites 0014 % involved in one compartment,then in the other compartment, non-proton 0015 % metabolites convert into non-proton reactants. In between is a reconstruction 0016 % transport reaction which must be elementally balanced for H to begin 0017 % with. We assume that the transporter involved has affinity only for a particular 0018 % species of metabolite defined by the reconstrution. 0019 % 0020 % INPUT 0021 % model.S 0022 % model.mets 0023 % model.SIntRxnBool Boolean of internal reactions 0024 % model.met(m).aveHbound average number of bound hydrogen ions 0025 % model.met(m).formula average number of bound hydrogen ions 0026 % model.met(m).charge charge of metabolite species from 0027 % reconstruction 0028 % 0029 % OPTIONAL INPUT 0030 % massImbalance nRxn x nElement matrix where massImbalance(i,j) is imbalance of 0031 % reaction i for element j. massImbalance(i,j)==0 if 0032 % balanced for that element. The first column is assumed 0033 % to correspond to H balancing. 0034 % 0035 % OUTPUT 0036 % model.S Stoichiometric matrix balanced for protons where each row. 0037 % corresponds to a reactant at specific pH. 0038 % 0039 % Ronan M. T. Fleming 0040 0041 if any(any(isnan(model.S))) 0042 error('NaN in S matrix before proton balancing.') 0043 end 0044 0045 [nMet,nRxn]=size(model.S); 0046 0047 %first save the reconstructions stoichiometric matrix 0048 model.Srecon = model.S; 0049 0050 if ~exist('massImbalance','var') 0051 massImBalancedBool=false(nRxn,1); 0052 else 0053 massImBalancedBool=(sum(massImbalance(:,2:end)')')~=0; 0054 end 0055 0056 %get the compartments of the model 0057 numChar=1; 0058 [allMetCompartments,uniqueMetCompartments]=getCompartment(model.mets,numChar); 0059 nUniqueCompartments=length(uniqueMetCompartments); 0060 0061 A=sparse(nMet,nUniqueCompartments); 0062 for m=1:nMet 0063 A(m,strcmp(allMetCompartments{m},uniqueMetCompartments))=1; 0064 end 0065 0066 compartmentHindex=zeros(length(uniqueMetCompartments),1); 0067 boolH=false(nMet,1); 0068 %indices of protons in different uniqueMetCompartments 0069 for p=1:length(uniqueMetCompartments) 0070 ind=find(strcmp(['h[' uniqueMetCompartments{p} ']'],model.mets)~=0); 0071 if ~isempty(ind) 0072 compartmentHindex(p)=ind; 0073 boolH(ind)=1; 0074 end 0075 end 0076 0077 %get the data on number of hydrogens bound by each reactant 0078 aveHbound=zeros(1,nMet); 0079 reconstructionH=zeros(1,nMet); 0080 metCharge=zeros(1,nMet); 0081 firstMissing=1; 0082 for m=1:nMet 0083 metCharge(m) = model.met(m).charge; 0084 aveHbound(1,m) = model.met(m).aveHbound; 0085 if strcmp(model.met(m).formula,'') 0086 if firstMissing 0087 warning('propper pH based balancing of protons requires a chemical formula for all metabolites') 0088 firstMissing=0; 0089 end 0090 reconstructionH(1,m)=NaN; 0091 else 0092 reconstructionH(1,m) = numAtomsOfElementInFormula(model.met(m).formula,'H'); 0093 end 0094 end 0095 0096 %check that canonical model is mass balanced 0097 dH = (reconstructionH*model.S)'; 0098 unbalancedInternalBool = dH~=0 & model.SIntRxnBool; 0099 0100 if any(unbalancedInternalBool) 0101 if 0 0102 fprintf('%s\n','Unbalanced reconstruction reactions:') 0103 unbalancedInd=find(unbalancedInternalBool); 0104 for p=1:length(unbalancedInd) 0105 fprintf('\n%20s\t%s\n',model.rxns{unbalancedInd(p)},model.rxn(unbalancedInd(p)).equation); 0106 end 0107 error('Hydrogen unbalanced reconstruction reactions exist.') 0108 else 0109 warning('Hydrogen unbalanced reconstruction reactions exist!') 0110 end 0111 end 0112 0113 %calculate the number of hydrogen atoms involved in each reaction 0114 Spos=model.S; 0115 Spos(model.S<0)=0; 0116 0117 %change in binding of Hydrogen ions which accompanies the conversion of 0118 %reactant into reconstruction metabolite species 0119 deltaHBound = aveHbound - reconstructionH; 0120 0121 fprintf('%s\n%s\n','Proton balancing of reactants.' ,'Assuming that transport reactions are specific to the metabolite species given in the reconstruction.') 0122 %assign new proton stoichiometric coefficients depending on compartment 0123 for n=1:nRxn 0124 if massImbalance(n,1)~=0 0125 warning('vonBertalanffy:pHbalanceProtons:OriginallyUnbalancedRxn', [model.rxns{n} ' reconstruction reaction not balanced for H to begin with']); % Changed from error to warning and added a message ID - Hulda 0126 end 0127 if strcmp('ACITL',model.rxns{n}) 0128 pause(eps) 0129 end 0130 0131 %no change for biomass reaction or exchange reactions or imbalanced 0132 %reactions 0133 if model.SIntRxnBool(n) && ~massImBalancedBool(n) 0134 %dont change reaction if any one of the metabolites has a NaN for 0135 %standard Gibbs energy of formation and therefore NaN for the 0136 %average number of H bound 0137 if ~any(isnan(aveHbound(model.S(:,n)~=0))) 0138 if 0 %debugging 0139 disp(model.rxns{n}) 0140 disp(model.rxn(n).equation) 0141 end 0142 %no change if only protons involved in a reaction 0143 if any(model.S(:,n)~=0 & ~boolH) 0144 0145 %uniqueMetCompartments of all metabolites involved in reactions 0146 metCompartments = allMetCompartments(model.S(:,n)~=0); 0147 rxnUniqueMetCompartments=unique(metCompartments); 0148 0149 %uniqueMetCompartments of non-proton metabolites involved in reactions 0150 metCompartmentsNotH = allMetCompartments(model.S(:,n)~=0 & ~boolH); 0151 uniqueCompartmentsNotH=unique(metCompartmentsNotH); 0152 0153 if length(rxnUniqueMetCompartments)>length(uniqueCompartmentsNotH) 0154 %proton transport across a membrane driving a reaction 0155 %make any necessary change to the proton stoichiometric 0156 %coefficient on the side with the rest of the 0157 %metabolites 0158 % e.g. abbreviation: 'ATPS4r' 0159 % officialName: 'ATP synthase (four protons for one ATP)' 0160 % equation: 'adp[c] + 4 h[e] + pi[c] <=> atp[c] + h2o[c] + 3 h[c] ' 0161 0162 %assumes the reconstruction transport reaction is elementally balanced for H 0163 compartBool = strcmp(uniqueCompartmentsNotH{1},allMetCompartments); 0164 metRxnCompartBool = model.S(:,n)~=0 & compartBool & ~boolH; 0165 % if 0 0166 % model.mets(model.S(:,n)~=0) 0167 % model.mets(metRxnCompartBool) 0168 % end 0169 % %index for H involved in compartment with reaction 0170 % indexHRxn = compartmentHindex(strcmp(uniqueCompartmentsNotH{1},uniqueMetCompartments)); 0171 0172 %find out if first index is substrate or product 0173 %compartment 0174 0175 % if sum(model.S(compartBool,n))<0 0176 % spIndexCol=1; 0177 % %if the hydrogen ion is a substrate, then store the number of 0178 % %hydrogen ions transported for the species reaction 0179 % %before mass balancing for the reactant reaction 0180 % netTransportZi(n,1)=-S(n,indexHRxn); 0181 % else 0182 % spIndexCol=2; 0183 % end 0184 % %save index of hydrogen ions for first compartment 0185 % substrateProductIndexH(n,spIndexCol)=indexHRxn; 0186 0187 %adjust the proton stoichiometric coefficient 0188 model.S(indexHRxn,n) = model.S(indexHRxn,n) - deltaHBound(metRxnCompartBool)*model.S(metRxnCompartBool,n); 0189 0190 % %index for H involved in compartment with reaction 0191 % indexHRxn = compartmentHindex(strcmp(uniqueCompartmentsNotH{1},uniqueMetCompartments)); 0192 % 0193 % %second index must be for opposite column of first 0194 % if spIndexCol==1 0195 % spIndexCol=2; 0196 % else 0197 % spIndexCol=1; 0198 % % %if the hydrogen ion is a substrate, then store the number of 0199 % % %hydrogen ions transported for the species reaction 0200 % % %before mass balancing for the reactant reaction 0201 % % netTransportZi(n,1)=-model.S(n,indexHRxn); 0202 % end 0203 % %save index of hydrogen ions for second compartment 0204 % substrateProductIndexH(n,spIndexCol)=indexHRxn; 0205 else 0206 %check the number of unique uniqueMetCompartments involving non 0207 %proton metabolites 0208 if length(uniqueCompartmentsNotH)==1 0209 %reaction involves metabolites in one compartment only 0210 indexHRxn = compartmentHindex(strcmp(uniqueCompartmentsNotH{1},uniqueMetCompartments)); 0211 %proton stoichiometric coefficient set to balance 0212 %protons with respect to average H bound by all other 0213 %substrates and products 0214 model.S(indexHRxn,n)= model.S(indexHRxn,n) - deltaHBound(~boolH)*model.S(~boolH,n); 0215 else 0216 %non-proton metabolites in two uniqueMetCompartments or more 0217 if length(uniqueCompartmentsNotH)>2 0218 error('More than two uniqueMetCompartments for a single reaction?!') 0219 end 0220 %model the transport across the membrane as three 0221 %reactions, one where non-proton reactants convert 0222 %into non-proton metabolites involved in one compartment, 0223 %then in the other compartment, non-proton 0224 %metabolites convert into non-proton reactants. In 0225 %between is a reconstruction transport reaction which must be 0226 %elementally balanced for H to begin with 0227 %this assumes that the transporter involved has 0228 %affinity only for a particular species of 0229 %metabolite defined by the reconstrution. 0230 0231 %first compartment 0232 compartBool1 = strcmp(uniqueCompartmentsNotH{1},allMetCompartments); 0233 %boolean for non-proton metabolites in first compartment 0234 metCompartBool1 = model.S(:,n)~=0 & compartBool1 & ~boolH; 0235 %index for stoichiometric coefficient of first compartment 0236 indexHRxn1 = compartmentHindex(strcmp(uniqueCompartmentsNotH{1},uniqueMetCompartments)); 0237 0238 % %find out if first index is substrate or product 0239 % %compartment 0240 % if sum(model.S(compartBool1,n))<0 0241 % spIndexCol=1; 0242 % else 0243 % spIndexCol=2; 0244 % end 0245 % %save index of hydrogen ion for first compartment 0246 % substrateProductIndexH(n,spIndexCol)=indexHRxn1; 0247 % 0248 % %net charge transported from first to second compartment 0249 % netTransportZi(n,1)=metCharge(compartBool1)*model.S(compartBool1,n); 0250 % 0251 % %TODO - need proper way to tell the order of uniqueMetCompartments 0252 % if model.S(indexHRxn1,n)~=0 0253 % if model.S(indexHRxn1,n)>0 0254 % warning('assuming it is a symport transport reaction - FIX') 0255 % netTransportZi(n,1)=-netTransportZi(n,1); 0256 % end 0257 % end 0258 0259 %second compartment 0260 compartBool2 = strcmp(uniqueCompartmentsNotH{2},allMetCompartments); 0261 %boolean for non-proton metabolites in first compartment 0262 metCompartBool2 = model.S(:,n)~=0 & compartBool2 & ~boolH; 0263 %index for stoichiometric coefficient of first compartment 0264 indexHRxn2 = compartmentHindex(strcmp(uniqueCompartmentsNotH{2},uniqueMetCompartments)); 0265 0266 % %second index must be for opposite column of first 0267 % if spIndexCol==1 0268 % spIndexCol=2; 0269 % else 0270 % spIndexCol=1; 0271 % end 0272 % %save index of hydrogen ion for second compartment 0273 % substrateProductIndexH(n,spIndexCol)=indexHRxn2; 0274 % 0275 % if 1 0276 % %net charge transported from second to first compartment 0277 % netTransportZiReverse=metCharge(compartBool2)*model.S(compartBool2,n); 0278 % if netTransportZiReverse~=netTransportZi(n,1) 0279 % error('Reconstruction reaction not charge balanced?'); 0280 % end 0281 % end 0282 0283 %mass balance reactant reactions 0284 model.S(indexHRxn1,n)=model.S(indexHRxn1,n) - deltaHBound(metCompartBool1)*model.S(metCompartBool1,n); 0285 model.S(indexHRxn2,n)=model.S(indexHRxn2,n) - deltaHBound(metCompartBool2)*model.S(metCompartBool2,n); 0286 0287 pause(eps) 0288 end 0289 end 0290 end 0291 end 0292 if aveHbound*model.S(:,n)>(eps*1e4) 0293 disp(model.rxns{n}) 0294 disp(model.rxn(n).equation) 0295 disp(aveHbound*model.S(:,n)) 0296 model.mets(model.S(:,n)~=0)' 0297 aveHbound(model.S(:,n)~=0) 0298 error(['Failure to proton balance. Reaction ' model.rxns{n} ', #' int2str(n)]) 0299 end 0300 if any(isnan(model.S(:,n))) 0301 pause(eps); 0302 end 0303 end 0304 if aveHbound*model.S(:,n)>(eps*1e4) 0305 if 1 %debugging 0306 disp(n) 0307 disp(model.rxns{n}) 0308 disp(model.rxn(n).equation) 0309 end 0310 error(['Reaction ' model.rxns{n} ', #' int2str(n) ' not proton balanced. No thermodynamic data available to balance reaction.']) 0311 end 0312 end 0313 if any(any(isnan(model.S))) 0314 error('NaN in S matrix after proton balancing.') 0315 end