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