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)
params parameter structure as used by getMILPParams (optional)
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,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 (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 % params parameter structure as used by getMILPParams (optional) 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 %If the user only supplied a single template model 0067 if ~iscell(models) 0068 models={models}; 0069 end 0070 0071 if nargin<3 0072 allowNetProduction=false; 0073 end 0074 if nargin<4 0075 useModelConstraints=false; 0076 end 0077 if nargin<5 0078 supressWarnings=false; 0079 end 0080 if nargin<6 0081 rxnScores=cell(numel(models),1); 0082 for i=1:numel(models) 0083 rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 0084 end 0085 end 0086 if isempty(rxnScores) 0087 rxnScores=cell(numel(models),1); 0088 for i=1:numel(models) 0089 rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 0090 end 0091 end 0092 if nargin<7 0093 params=struct(); 0094 end 0095 0096 if ~iscell(rxnScores) 0097 rxnScores={rxnScores}; 0098 end 0099 0100 models=models(:); 0101 rxnScores=rxnScores(:); 0102 0103 %Check if the original model has an unconstrained field. If so, give a 0104 %warning 0105 if supressWarnings==false 0106 if isfield(model,'unconstrained') 0107 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'; 0108 dispEM(EM,false); 0109 else 0110 if isempty(getExchangeRxns(model,'both')) 0111 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'); 0112 end 0113 end 0114 end 0115 0116 %Simplify the template models to remove constrained rxns. At the same time, 0117 %check that the id of the template models isn't the same as the model. That 0118 %would cause an error further down 0119 for i=1:numel(models) 0120 models{i}.rxnScores=rxnScores{i}; 0121 models{i}=simplifyModel(models{i},false,false,true); 0122 if strcmpi(models{i}.id,model.id) 0123 EM='The reference model(s) cannot have the same id as the model'; 0124 dispEM(EM); 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 %removeReactions 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],'metNames',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 exitFlag=1; 0154 return; 0155 end 0156 0157 %Then check that the merged model is solveable 0158 sol=solveLP(allModels); 0159 if isempty(sol.f) 0160 EM='There are no reactions in the template model(s) that can make the model constraints satisfied'; 0161 dispEM(EM); 0162 end 0163 0164 %Remove dead ends for speed reasons. This has to be done here and 0165 %duplicate below because there is otherwise a risk that a reaction 0166 %which is constrained to something relevant is removed 0167 allModels=simplifyModel(allModels,false,false,false,true,false,false,false,[],true); 0168 allModels.c(:)=0; 0169 else 0170 %Remove dead ends for speed reasons 0171 allModels=simplifyModel(allModels,false,false,false,true,false,false,false,[],true); 0172 allModels.c(:)=0; 0173 0174 %If model constraints shouldn't be used, then determine which reactions 0175 %to force to have flux 0176 0177 %Get the reactions that can carry flux in the original model 0178 originalFlux=haveFlux(model,1); 0179 0180 %For the ones that can't carry flux, see if they can do so in the 0181 %merged model 0182 toCheck=intersect(allModels.rxns(strcmp(allModels.rxnFrom,model.id)),model.rxns(~originalFlux)); 0183 0184 %Get the ones that still cannot carry flux 0185 I=haveFlux(allModels,1,toCheck); 0186 0187 %Get the reactions that can't carry flux in the original model, but can 0188 %in the merged one 0189 K=toCheck(I); 0190 0191 %This is a temporary thing to only look at the non-reversible rxns. 0192 %This is because all reversible rxns can have a flux in the 0193 %irreversible model format that is used by getMinNrFluxes 0194 [~, I]=ismember(K,model.rxns); 0195 K(model.rev(I)~=0)=[]; 0196 0197 %Constrain all reactions in the original model to have a flux 0198 allModels.lb(ismember(allModels.rxns,K))=0.1; 0199 0200 %Return stuff 0201 newConnected=K; 0202 cannotConnect=setdiff(model.rxns(~originalFlux),newConnected); 0203 end 0204 0205 %Then minimize for the number of fluxes used. The fixed rxns doesn't need 0206 %to participate 0207 templateRxns=find(~strcmp(allModels.rxnFrom,model.id)); 0208 [~, J, exitFlag]=getMinNrFluxes(allModels,templateRxns,params,allModels.rxnScores(templateRxns)); 0209 0210 %Remove everything except for the added ones 0211 I=true(numel(allModels.rxns),1); 0212 I(templateRxns(J))=false; 0213 addedModel=removeReactions(allModels,I,true,true,true); 0214 0215 newModel=mergeModels({model;addedModel},'metNames',true); 0216 addedRxns=setdiff(newModel.rxns,model.rxns); 0217 newModel=rmfield(newModel,'rxnScores'); 0218 end