runINIT
Generates a model using the INIT algorithm, based on proteomics and/or
transcriptomics and/or metabolomics and/or metabolic tasks. This is the
original implementation, which is now replaced with ftINIT.
model a reference model structure
rxnScores a vector of scores for the reactions in the model.
Positive scores are reactions to keep and negative
scores are reactions to exclude (optional, default all 0.0)
presentMets cell array with unique metabolite names that the model
should produce (optional, default [])
essentialRxns cell array of reactions that are essential and that
have to be in the resulting model. This is normally
used when fitting a model to task (see fitTasks) (optional,
default [])
prodWeight a score that determines the value of having
net-production of metabolites. This is a way of having
a more functional network as it provides a reason for
including bad reactions for connectivity reasons. This
score is for each metabolite, and the sum of these weights
and the scores for the reactions is what is optimized
(optional, default 0.5)
allowExcretion true if excretion of all metabolites should be allowed.
This results in fewer reactions being considered
dead-ends, but all reactions in the resulting model may
not be able to carry flux. If this is "false" then the
equality constraints are taken from model.b. If the
input model lacks exchange reactions then this should
probably be "true", or a large proportion of the model
would be excluded for connectivity reasons
(optional, default false)
noRevLoops true if reversible reactions should be constrained to
only carry flux in one direction. This prevents
reversible reactions from being wrongly assigned as
connected (the forward and backward reactions can form a
loop and therefore appear connected), but it makes the
problem significantly more computationally intensive to
solve (two more integer constraints per reversible reaction)
(optional, default false)
outModel the resulting model structure
deletedRxns reactions which were deleted by the algorithm
metProduction array that indicates which of the
metabolites in presentMets that could be
produced
-2: metabolite name not found in model
-1: metabolite found, but it could not be produced
1: metabolite could be produced
fValue objective value (sum of (the negative of)
reaction scores for the included reactions and
prodWeight*number of produced metabolites)
This function is the actual implementation of the algorithm. See
getINITModel for a higher-level function for model reconstruction. See
PLoS Comput Biol. 2012;8(5):e1002518 for details regarding the
implementation.
Usage: [outModel deletedRxns metProduction fValue]=runINIT(model,...
rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,...
noRevLoops,params)0001 function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops) 0002 % runINIT 0003 % Generates a model using the INIT algorithm, based on proteomics and/or 0004 % transcriptomics and/or metabolomics and/or metabolic tasks. This is the 0005 % original implementation, which is now replaced with ftINIT. 0006 % 0007 % model a reference model structure 0008 % rxnScores a vector of scores for the reactions in the model. 0009 % Positive scores are reactions to keep and negative 0010 % scores are reactions to exclude (optional, default all 0.0) 0011 % presentMets cell array with unique metabolite names that the model 0012 % should produce (optional, default []) 0013 % essentialRxns cell array of reactions that are essential and that 0014 % have to be in the resulting model. This is normally 0015 % used when fitting a model to task (see fitTasks) (optional, 0016 % default []) 0017 % prodWeight a score that determines the value of having 0018 % net-production of metabolites. This is a way of having 0019 % a more functional network as it provides a reason for 0020 % including bad reactions for connectivity reasons. This 0021 % score is for each metabolite, and the sum of these weights 0022 % and the scores for the reactions is what is optimized 0023 % (optional, default 0.5) 0024 % allowExcretion true if excretion of all metabolites should be allowed. 0025 % This results in fewer reactions being considered 0026 % dead-ends, but all reactions in the resulting model may 0027 % not be able to carry flux. If this is "false" then the 0028 % equality constraints are taken from model.b. If the 0029 % input model lacks exchange reactions then this should 0030 % probably be "true", or a large proportion of the model 0031 % would be excluded for connectivity reasons 0032 % (optional, default false) 0033 % noRevLoops true if reversible reactions should be constrained to 0034 % only carry flux in one direction. This prevents 0035 % reversible reactions from being wrongly assigned as 0036 % connected (the forward and backward reactions can form a 0037 % loop and therefore appear connected), but it makes the 0038 % problem significantly more computationally intensive to 0039 % solve (two more integer constraints per reversible reaction) 0040 % (optional, default false) 0041 % 0042 % outModel the resulting model structure 0043 % deletedRxns reactions which were deleted by the algorithm 0044 % metProduction array that indicates which of the 0045 % metabolites in presentMets that could be 0046 % produced 0047 % -2: metabolite name not found in model 0048 % -1: metabolite found, but it could not be produced 0049 % 1: metabolite could be produced 0050 % fValue objective value (sum of (the negative of) 0051 % reaction scores for the included reactions and 0052 % prodWeight*number of produced metabolites) 0053 % 0054 % This function is the actual implementation of the algorithm. See 0055 % getINITModel for a higher-level function for model reconstruction. See 0056 % PLoS Comput Biol. 2012;8(5):e1002518 for details regarding the 0057 % implementation. 0058 % 0059 % Usage: [outModel deletedRxns metProduction fValue]=runINIT(model,... 0060 % rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,... 0061 % noRevLoops,params) 0062 0063 if nargin<2 0064 rxnScores=zeros(numel(model.rxns),1); 0065 end 0066 if isempty(rxnScores) 0067 rxnScores=zeros(numel(model.rxns),1); 0068 end 0069 if nargin<3 || isempty(presentMets) 0070 presentMets={}; 0071 else 0072 presentMets=convertCharArray(presentMets); 0073 end 0074 if nargin<4 || isempty(essentialRxns) 0075 essentialRxns={}; 0076 else 0077 essentialRxns=convertCharArray(essentialRxns); 0078 end 0079 if nargin<5 || isempty(prodWeight) 0080 prodWeight=0.5; 0081 end 0082 if nargin<6 0083 allowExcretion=false; 0084 end 0085 if nargin<7 0086 noRevLoops=false; 0087 end 0088 if nargin<8 0089 params=[]; 0090 end 0091 0092 if numel(presentMets)~=numel(unique(presentMets)) 0093 EM='Duplicate metabolite names in presentMets'; 0094 dispEM(EM); 0095 end 0096 0097 %Default is that the metabolites cannot be produced 0098 if ~isempty(presentMets) 0099 metProduction=ones(numel(presentMets),1)*-2; 0100 presentMets=upper(presentMets); 0101 pmIndexes=find(ismember(presentMets,upper(model.metNames))); 0102 metProduction(pmIndexes)=-1; %Then set that they are at least found 0103 else 0104 metProduction=[]; 0105 pmIndexes=[]; 0106 end 0107 0108 %The model should be in the reversible format and all relevant exchange 0109 %reactions should be open 0110 if isfield(model,'unconstrained') 0111 EM='Exchange metabolites are still present in the model. Use simplifyModel if this is not intended'; 0112 dispEM(EM,false); 0113 end 0114 0115 %The irreversible reactions that are essential must have a flux and are 0116 %therefore not optimized for using MILP, which reduces the problem size. 0117 %However, reversible reactions must have a flux in one direction, so they 0118 %have to stay in the problem. The essentiality constraint on reversible 0119 %reactions is implemented in the same manner as for reversible reactions 0120 %when noRevLoops==true, but with the additional constraint that C ub=-1. 0121 %This forces one of the directions to be active. 0122 revRxns=find(model.rev~=0); 0123 essentialReversible=find(ismember(model.rxns(revRxns),essentialRxns)); 0124 essentialRxns=intersect(essentialRxns,model.rxns(model.rev==0)); 0125 0126 %Convert the model to irreversible 0127 irrevModel=convertToIrrev(model); 0128 rxnScores=[rxnScores;rxnScores(model.rev==1)]; 0129 %These are used if noRevLoops is true 0130 if noRevLoops==true 0131 forwardIndexes=find(model.rev~=0); 0132 backwardIndexes=(numel(model.rxns)+1:numel(irrevModel.rxns))'; 0133 else 0134 %Then they should only be used for essential reversible reactions 0135 forwardIndexes=revRxns(essentialReversible); 0136 backwardIndexes=essentialReversible+numel(model.rxns); 0137 end 0138 0139 %Get the indexes of the essential reactions and remove them from the 0140 %scoring vector 0141 essentialIndex=find(ismember(irrevModel.rxns,essentialRxns)); 0142 rxnScores(essentialIndex)=[]; 0143 0144 %Go through each of the presentMets (if they exist) and modify the S matrix 0145 %so that each reaction which produces any of them also produces a 0146 %corresponding fake metabolite and the opposite in the reverse direction. 0147 0148 %This is to deal with the fact that there is no compartment info regarding 0149 %the presentMets. This modifies the irrevModel structure, but that is fine 0150 %since it's the model structure that is returned. 0151 if any(pmIndexes) 0152 irrevModel.metNames=upper(irrevModel.metNames); 0153 metsToAdd.mets=strcat({'FAKEFORPM'},num2str(pmIndexes)); 0154 metsToAdd.metNames=metsToAdd.mets; 0155 metsToAdd.compartments=irrevModel.comps{1}; 0156 0157 %There is no constraints on the metabolites yet, since maybe not all of 0158 %them could be produced 0159 irrevModel=addMets(irrevModel,metsToAdd); 0160 end 0161 0162 %Modify the matrix 0163 for i=1:numel(pmIndexes) 0164 %Get the matching mets 0165 I=ismember(irrevModel.metNames,presentMets(pmIndexes(i))); 0166 0167 %Find the reactions where any of them are used. 0168 [~, K, L]=find(irrevModel.S(I,:)); 0169 0170 %This ugly loop is to avoid problems if a metabolite occurs several 0171 %times in one reaction 0172 KK=unique(K); 0173 LL=zeros(numel(KK),1); 0174 for j=1:numel(KK) 0175 LL(j)=sum(L(K==KK(j))); 0176 end 0177 irrevModel.S(numel(irrevModel.mets)-numel(pmIndexes)+i,KK)=LL; 0178 end 0179 0180 %Some nice to have numbers 0181 nMets=numel(irrevModel.mets); 0182 nRxns=numel(irrevModel.rxns); 0183 nEssential=numel(essentialIndex); 0184 nNonEssential=nRxns-nEssential; 0185 nonEssentialIndex=setdiff(1:nRxns,essentialIndex); 0186 S=irrevModel.S; 0187 0188 %Add so that each non-essential reaction produces one unit of a fake 0189 %metabolite 0190 temp=sparse(1:nRxns,1:nRxns,1); 0191 temp(essentialIndex,:)=[]; 0192 S=[S;temp]; 0193 0194 %Add another set of reactions (will be binary) which also produce these 0195 %fake metabolites, but with a stoichiometry of 1000 0196 temp=sparse(1:nNonEssential,1:nNonEssential,1000); 0197 temp=[sparse(nMets,nNonEssential);temp]; 0198 S=[S temp]; 0199 0200 %Add reactions for net-production of (real) metabolites 0201 if prodWeight~=0 0202 temp=[speye(nMets-numel(pmIndexes))*-1;sparse(nNonEssential+numel(pmIndexes),nMets-numel(pmIndexes))]; 0203 S=[S temp]; 0204 %To keep the number of reactions added like this 0205 nNetProd=nMets-numel(pmIndexes); 0206 else 0207 nNetProd=0; 0208 end 0209 0210 %Add constraints so that reversible reactions can only be used in one 0211 %direction. This is done by adding the fake metabolites A, B, C for each 0212 %reversible reaction in the following manner 0213 % forward: A + .. => ... backwards: B + ... => ... int1: C => 1000 A int2: 0214 % C => 1000 B A ub=999.9 B ub=999.9 C lb=-1 int1 and int2 are binary 0215 if any(forwardIndexes) 0216 nRevBounds=numel(forwardIndexes); 0217 0218 %Add the A metabolites for the forward reactions and the B metabolites 0219 %for the reverse reactions 0220 I=speye(numel(irrevModel.rxns))*-1; 0221 temp=[I(forwardIndexes,:);I(backwardIndexes,:)]; 0222 0223 %Padding 0224 temp=[temp sparse(size(temp,1),size(S,2)-numel(irrevModel.rxns))]; 0225 0226 %Add the int1 & int2 reactions that produce A and B 0227 temp=[temp speye(nRevBounds*2)*1000]; 0228 0229 %And add that they also consume C 0230 temp=[temp;[sparse(nRevBounds,size(S,2)) speye(nRevBounds)*-1 speye(nRevBounds)*-1]]; 0231 0232 %Add the new reactions and metabolites 0233 S=[S sparse(size(S,1),nRevBounds*2)]; 0234 S=[S;temp]; 0235 else 0236 nRevBounds=0; 0237 end 0238 0239 %Add so that the essential reactions must have a small flux and that the 0240 %binary ones (and net-production reactions) may have zero flux. The integer 0241 %reactions for reversible reactions have [0 1] 0242 prob.blx=[irrevModel.lb;zeros(nNonEssential+nNetProd+nRevBounds*2,1)]; 0243 prob.blx(essentialIndex)=max(0.1,prob.blx(essentialIndex)); 0244 0245 %Add so that the binary ones and net-production reactions can have at the 0246 %most flux 1.0 0247 prob.bux=[irrevModel.ub;ones(nNonEssential+nNetProd+nRevBounds*2,1)]; 0248 0249 %Add that the fake metabolites must be produced in a small amount and that 0250 %the A and B metabolites for reversible reactions can be [0 999.9] and C 0251 %metabolites [-1 0] 0252 prob.blc=[irrevModel.b(:,1);ones(nNonEssential,1);zeros(nRevBounds*2,1);ones(nRevBounds,1)*-1]; 0253 0254 %Add that normal metabolites can be freely excreted if 0255 %allowExcretion==true, and that the fake ones can be excreted 1000 units at 0256 %most. C metabolites for essential reversible reactions should have an 0257 %upper bound of -1. If noRevLoops is false, then add this constraint for 0258 %all the reactions instead. 0259 if noRevLoops==true 0260 revUB=zeros(nRevBounds,1); 0261 revUB(essentialReversible)=-1; 0262 else 0263 revUB=ones(nRevBounds,1)*-1; 0264 end 0265 if allowExcretion==true 0266 metUB=inf(nMets,1); 0267 else 0268 metUB=irrevModel.b(:,min(size(irrevModel.b,2),2)); 0269 end 0270 prob.buc=[metUB;ones(nNonEssential,1)*1000;ones(nRevBounds*2,1)*999.9;revUB]; 0271 0272 %Add objective coefficients for the binary reactions. The negative is used 0273 %since we're minimizing. The negative is taken for the prodWeight as well, 0274 %in order to be consistent with the syntax that positive scores are good 0275 prob.c=[zeros(nRxns,1);rxnScores;ones(nNetProd,1)*prodWeight*-1;zeros(nRevBounds*2,1)]; 0276 prob.a=S; 0277 0278 % adapt problem structure for cobra-style solver 0279 prob.c=[prob.c;zeros(size(prob.a,1),1)]; 0280 prob.A=[prob.a -speye(size(prob.a,1))]; 0281 prob.b=zeros(size(prob.a,1), 1); 0282 prob.ub=[prob.bux; prob.buc]; 0283 prob.osense=1; 0284 prob.csense=char(zeros(1,size(prob.a,1))); 0285 prob.csense(:)='E'; 0286 0287 %We still don't know which of the presentMets that can be produced. Go 0288 %through them, force production, and see if the problem can be solved 0289 for i=1:numel(pmIndexes) 0290 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=1; 0291 prob.lb=[prob.blx; prob.blc]; 0292 res=optimizeProb(prob,params); 0293 isFeasible=checkSolution(res); 0294 if ~isFeasible 0295 %Reset the constraint again 0296 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=0; 0297 else 0298 %Metabolite produced 0299 metProduction(pmIndexes(i))=1; 0300 end 0301 end 0302 prob.lb=[prob.blx; prob.blc]; 0303 0304 %Add that the binary reactions may only take integer values. 0305 prob.vartype = repmat('C', 1, size(prob.A, 2)); 0306 allInt=[(nRxns+1):(nRxns+nNonEssential) size(S,2)-nRevBounds*2+1:size(S,2)]; 0307 prob.vartype(allInt) = 'B'; 0308 0309 % solve problem 0310 res=optimizeProb(prob,params); 0311 0312 %Problem should not be infeasible, but it is possible that the time limit 0313 %was reached before finding any solutions. 0314 if ~checkSolution(res) 0315 if strcmp(res.origStat, 'TIME_LIMIT') 0316 EM='Time limit reached without finding a solution. Try increasing the TimeLimit parameter.'; 0317 else 0318 EM='The problem is infeasible'; 0319 end 0320 dispEM(EM); 0321 end 0322 0323 fValue=res.obj; 0324 0325 %Get all reactions used in the irreversible model 0326 usedRxns=(nonEssentialIndex(res.full(nRxns+1:nRxns+nNonEssential)<0.1))'; 0327 0328 %Map to reversible model IDs 0329 usedRxns=[usedRxns(usedRxns<=numel(model.rxns));revRxns(usedRxns(usedRxns>numel(model.rxns))-numel(model.rxns))]; 0330 0331 %Then get the ones that are not used in either direction or is essential 0332 I=true(numel(model.rxns),1); 0333 I(usedRxns)=false; 0334 I(essentialIndex)=false; 0335 deletedRxns=model.rxns(I); 0336 outModel=removeReactions(model,I,true,true); 0337 end