fillGaps Uses template model(s) to fill gaps in a model model a model structure that may contains gaps to be filled models a cell array of reference models or a model structure. The gaps will be filled using reactions from these models allowNetProduction true if net production of all metabolites is allowed. A reaction can be unable to carry flux because one of the reactants is unavailable or because one of the products can't be further processed. If this parameter is true, only the first type of unconnectivity is considered (opt, default false) useModelConstraints true if the constraints specified in the model structure should be used. If false then reactions included from the template model(s) so that as many reactions as possible in model can carry flux (opt, default false) supressWarnings false if warnings should be displayed (opt, default false) rxnScores array with scores for each of the reactions in the reference model(s). If more than one model is supplied, then rxnScores should be a cell array of vectors. The solver will try to maximize the sum of the scores for the included reactions (opt, default is -1 for all reactions) params parameter structure as used by getMILPParams (opt) newConnected cell array with the reactions that could be connected. This is not calulated if useModelConstraints is true cannotConnect cell array with reactions that could not be connected. This is not calculated if useModelConstraints is true addedRxns cell array with the reactions that were added from "models" newModel the model with reactions added to fill gaps exitFlag 1: optimal solution found -1: no feasible solution found -2: optimization time out This method works by merging the model to the reference model(s) and checking which reactions can carry flux. All reactions that can't carry flux are removed (cannotConnect). If useModelConstraints is false it then solves the MILP problem of minimizing the number of active reactions from the reference models that are required to have flux in all the reactions in model. This requires that the input model has exchange reactions present for the nutrients that are needed for its metabolism. If useModelConstraints is true then the problem is to include as few reactions as possible from the reference models in order to satisfy the model constraints. The intended use is that the user can attempt a general gap-filling using useModelConstraint=false or a more targeted gap-filling by setting constraints in the model structure and then use useModelConstraints=true. Say that the user want to include reactions so that all biomass components can be synthesized. He/she could then define a biomass equation and set the lower bound to >0. Running this function with useModelConstraints=true would then give the smallest set of reactions that have to be included in order for the model to produce biomass. Usage: [newConnected cannotConnect addedRxns newModel exitFlag]=... fillGaps(model,models,allowNetProduction,useModelConstraints,... supressWarnings,rxnScores,params) Rasmus Agren, 2013-07-16
0001 function [newConnected cannotConnect addedRxns newModel exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params) 0002 % fillGaps 0003 % Uses template model(s) to fill gaps in a model 0004 % 0005 % model a model structure that may contains gaps to be filled 0006 % models a cell array of reference models or a model structure. 0007 % The gaps will be filled using reactions from these models 0008 % allowNetProduction true if net production of all metabolites is 0009 % allowed. A reaction can be unable to carry flux because one of 0010 % the reactants is unavailable or because one of the 0011 % products can't be further processed. If this 0012 % parameter is true, only the first type of 0013 % unconnectivity is considered (opt, default false) 0014 % useModelConstraints true if the constraints specified in the model 0015 % structure should be used. If false then reactions 0016 % included from the template model(s) so that as many 0017 % reactions as possible in model can carry flux 0018 % (opt, default false) 0019 % supressWarnings false if warnings should be displayed (opt, default 0020 % false) 0021 % rxnScores array with scores for each of the reactions in the 0022 % reference model(s). If more than one model is supplied, 0023 % then rxnScores should be a cell array of vectors. 0024 % The solver will try to maximize the sum of the 0025 % scores for the included reactions (opt, default 0026 % is -1 for all reactions) 0027 % params parameter structure as used by getMILPParams (opt) 0028 % 0029 % newConnected cell array with the reactions that could be 0030 % connected. This is not calulated if 0031 % useModelConstraints is true 0032 % cannotConnect cell array with reactions that could not be 0033 % connected. This is not calculated if 0034 % useModelConstraints is true 0035 % addedRxns cell array with the reactions that were added from 0036 % "models" 0037 % newModel the model with reactions added to fill gaps 0038 % exitFlag 1: optimal solution found 0039 % -1: no feasible solution found 0040 % -2: optimization time out 0041 % 0042 % This method works by merging the model to the reference model(s) and 0043 % checking which reactions can carry flux. All reactions that can't 0044 % carry flux are removed (cannotConnect). 0045 % If useModelConstraints is false it then solves the MILP problem of 0046 % minimizing the number of active reactions from the reference models 0047 % that are required to have flux in all the reactions in model. This 0048 % requires that the input model has exchange reactions present for the 0049 % nutrients that are needed for its metabolism. If useModelConstraints is 0050 % true then the problem is to include as few reactions as possible from 0051 % the reference models in order to satisfy the model constraints. 0052 % The intended use is that the user can attempt a general gap-filling using 0053 % useModelConstraint=false or a more targeted gap-filling by setting 0054 % constraints in the model structure and then use 0055 % useModelConstraints=true. Say that the user want to include reactions 0056 % so that all biomass components can be synthesized. He/she could then 0057 % define a biomass equation and set the lower bound to >0. Running this 0058 % function with useModelConstraints=true would then give the smallest set 0059 % of reactions that have to be included in order for the model to produce 0060 % biomass. 0061 % 0062 % Usage: [newConnected cannotConnect addedRxns newModel exitFlag]=... 0063 % fillGaps(model,models,allowNetProduction,useModelConstraints,... 0064 % supressWarnings,rxnScores,params) 0065 % 0066 % Rasmus Agren, 2013-07-16 0067 % 0068 0069 %If the user only supplied a single template model 0070 if ~iscell(models) 0071 models={models}; 0072 end 0073 0074 if nargin<3 0075 allowNetProduction=false; 0076 end 0077 if nargin<4 0078 useModelConstraints=false; 0079 end 0080 if nargin<5 0081 supressWarnings=false; 0082 end 0083 if nargin<6 0084 rxnScores=cell(numel(models),1); 0085 for i=1:numel(models) 0086 rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 0087 end 0088 end 0089 if isempty(rxnScores) 0090 rxnScores=cell(numel(models),1); 0091 for i=1:numel(models) 0092 rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 0093 end 0094 end 0095 if nargin<7 0096 params=[]; 0097 end 0098 0099 if ~iscell(rxnScores) 0100 rxnScores={rxnScores}; 0101 end 0102 0103 models=models(:); 0104 rxnScores=rxnScores(:); 0105 0106 %Check if the original model has an unconstrained field. If so, give a warning 0107 if supressWarnings==false 0108 if isfield(model,'unconstrained'); 0109 fprintf('WARNING: This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model still has the "unconstrained" field.\n'); 0110 else 0111 if isempty(getExchangeRxns(model,'both')) 0112 fprintf('NOTE: This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model does not seem to contain any such reactions.\n'); 0113 end 0114 end 0115 end 0116 0117 %Simplify the template models to remove constrained rxns. At the same time, 0118 %check that the id of the template models isn't the same as the model. 0119 %That would cause an error further down 0120 for i=1:numel(models) 0121 models{i}.rxnScores=rxnScores{i}; 0122 models{i}=simplifyModel(models{i},false,false,true); 0123 if strcmpi(models{i}.id,model.id) 0124 throw(MException('','The reference model(s) cannot have the same id as the model')); 0125 end 0126 end 0127 0128 %This is a rather ugly solution to the issue that it's a bit tricky to keep 0129 %track of which scores belong to which reactions. This requires that 0130 %removeRxns and mergeModels are modified to check for the new field. 0131 model.rxnScores=zeros(numel(model.rxns),1); 0132 0133 %First merge all models into one big one 0134 allModels=mergeModels([{model};models],true); 0135 0136 %Add that net production is ok 0137 if allowNetProduction==true 0138 %A second column in model.b means that the b field is lower and upper 0139 %bound on the RHS. 0140 model.b=[model.b(:,1) inf(numel(model.mets),1)]; 0141 allModels.b=[allModels.b(:,1) inf(numel(allModels.mets),1)]; 0142 end 0143 0144 if useModelConstraints==true 0145 newConnected={}; 0146 cannotConnect={}; 0147 0148 %Check that the input model isn't solveable without any input 0149 sol=solveLP(model); 0150 if ~isempty(sol.f) 0151 addedRxns={}; 0152 newModel=model; 0153 return; 0154 end 0155 0156 %Then check that the merged model is solveable 0157 sol=solveLP(allModels); 0158 if isempty(sol.f) 0159 throw(MException('','There are no reactions in the template model(s) that can make the model constraints satisfied')); 0160 end 0161 0162 %Remove dead ends for speed reasons. This has to be done here and 0163 %duplicate below because there is otherwise a risk that a reaction 0164 %which is constrained to something relevant is removed 0165 allModels=simplifyModel(allModels,false,false,false,true,false,false,[],true); 0166 allModels.c(:)=0; 0167 else 0168 %Remove dead ends for speed reasons 0169 allModels=simplifyModel(allModels,false,false,false,true,false,false,[],true); 0170 allModels.c(:)=0; 0171 0172 %If model constraints shouldn't be used, then determine which reactions to 0173 %force to have flux 0174 0175 %Get the reactions that can carry flux in the original model 0176 originalFlux=haveFlux(model,1); 0177 0178 %For the ones that can't carry flux, see if they can do so in the merged 0179 %model 0180 toCheck=intersect(allModels.rxns(strcmp(allModels.rxnFrom,model.id)),model.rxns(~originalFlux)); 0181 0182 %Get the ones that still cannot carry flux 0183 I=haveFlux(allModels,1,toCheck); 0184 0185 %Get the reactions that can't carry flux in the original model, but can 0186 %in the merged one 0187 K=toCheck(I); 0188 0189 %This is a temporary thing to only look at the non-reversible rxns. 0190 %This is because all reversible rxns can have a flux in the irreversible 0191 %model format that is used by getMinNrFluxes 0192 [crap I]=ismember(K,model.rxns); 0193 K(model.rev(I)~=0)=[]; 0194 0195 %Constrain all reactions in the original model to have a flux 0196 allModels.lb(ismember(allModels.rxns,K))=0.1; 0197 0198 %Return stuff 0199 newConnected=K; 0200 cannotConnect=setdiff(model.rxns(~originalFlux ),newConnected); 0201 end 0202 0203 %Then minimize for the number of fluxes used. The fixed rxns doesn't need 0204 %to participate 0205 templateRxns=find(~strcmp(allModels.rxnFrom,model.id)); 0206 [crap J exitFlag]=getMinNrFluxes(allModels,templateRxns,params,allModels.rxnScores(templateRxns)); 0207 0208 %Remove everything except for the added ones 0209 I=true(numel(allModels.rxns),1); 0210 I(templateRxns(J))=false; 0211 addedModel=removeRxns(allModels,I,true); 0212 0213 newModel=mergeModels({model;addedModel},true); 0214 addedRxns=setdiff(newModel.rxns,model.rxns); 0215 newModel=rmfield(newModel,'rxnScores'); 0216 end