pHbalanceProtons

PURPOSE ^

Mass balance protons for each reaction by adjusting the hydrogen ion stoichiometric coefficient.

SYNOPSIS ^

function model=pHbalanceProtons(model,massImbalance)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003