ftINITFillGaps Variant of fillGaps specially adapted to speed up generation of ftINIT models. tModel model that contains the task-specific rxns. origModel model without task-specific rxns. tRefModel reference tModel - the full tModel, containing all rxns used for gap-filling of tModel + the task-specific rxns 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 tModel structure should be used. If false then reactions included from the template tModel(s) so that as many reactions as possible in tModel can carry flux (optional, default false) supressWarnings false if warnings should be displayed (optional, default false) rxnScores scores for each of the reactions in the reference tModel. The solver will try to maximize the sum of the scores for the included reactions params parameter structure as used by getMILPParams verbose if true, the MILP progression will be shown. addedRxns the rxns added newModel the tModel 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 tModel to the reference model and checking which reactions can carry flux. All reactions that can't carry flux are removed. If useModelConstraints is false it then solves the MILP problem of minimizing the number of active reactions from the reference model that are required to have flux in all the reactions in model. This requires that the input tModel 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 tModel 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 tModel to produce biomass. Usage: [newModel, exitFlag]=... fillGaps(tModel,models,allowNetProduction,... supressWarnings,rxnScores,params)
0001 function [addedRxns, newModel, exitFlag]=ftINITFillGaps(tModel, origModel, tRefModel,allowNetProduction,supressWarnings,rxnScores,params,verbose) 0002 % ftINITFillGaps 0003 % Variant of fillGaps specially adapted to speed up generation of ftINIT models. 0004 % 0005 % tModel model that contains the task-specific rxns. 0006 % origModel model without task-specific rxns. 0007 % tRefModel reference tModel - the full tModel, containing all rxns 0008 % used for gap-filling of tModel + the task-specific rxns 0009 % allowNetProduction true if net production of all metabolites is 0010 % allowed. A reaction can be unable to carry flux because one of 0011 % the reactants is unavailable or because one of the 0012 % products can't be further processed. If this 0013 % parameter is true, only the first type of 0014 % unconnectivity is considered (optional, default false) 0015 % useModelConstraints true if the constraints specified in the tModel 0016 % structure should be used. If false then reactions 0017 % included from the template tModel(s) so that as many 0018 % reactions as possible in tModel can carry flux 0019 % (optional, default false) 0020 % supressWarnings false if warnings should be displayed (optional, default 0021 % false) 0022 % rxnScores scores for each of the reactions in the 0023 % reference tModel. 0024 % The solver will try to maximize the sum of the 0025 % scores for the included reactions 0026 % params parameter structure as used by getMILPParams 0027 % verbose if true, the MILP progression will be shown. 0028 % 0029 % addedRxns the rxns added 0030 % newModel the tModel with reactions added to fill gaps 0031 % exitFlag 1: optimal solution found 0032 % -1: no feasible solution found 0033 % -2: optimization time out 0034 % 0035 % This method works by merging the tModel to the reference model and 0036 % checking which reactions can carry flux. All reactions that can't 0037 % carry flux are removed. 0038 % If useModelConstraints is false it then solves the MILP problem of 0039 % minimizing the number of active reactions from the reference model 0040 % that are required to have flux in all the reactions in model. This 0041 % requires that the input tModel has exchange reactions present for the 0042 % nutrients that are needed for its metabolism. If useModelConstraints is 0043 % true then the problem is to include as few reactions as possible from 0044 % the reference models in order to satisfy the tModel constraints. 0045 % The intended use is that the user can attempt a general gap-filling using 0046 % useModelConstraint=false or a more targeted gap-filling by setting 0047 % constraints in the model structure and then use 0048 % useModelConstraints=true. Say that the user want to include reactions 0049 % so that all biomass components can be synthesized. He/she could then 0050 % define a biomass equation and set the lower bound to >0. Running this 0051 % function with useModelConstraints=true would then give the smallest set 0052 % of reactions that have to be included in order for the tModel to produce 0053 % biomass. 0054 % 0055 % Usage: [newModel, exitFlag]=... 0056 % fillGaps(tModel,models,allowNetProduction,... 0057 % supressWarnings,rxnScores,params) 0058 0059 if isempty(rxnScores) 0060 rxnScores=cell(numel(models),1); 0061 for i=1:numel(models) 0062 rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 0063 end 0064 end 0065 0066 %Simplify the template models to remove constrained rxns. At the same time, 0067 %check that the id of the template models isn't the same as the tModel. That 0068 %would cause an error further down 0069 tRefModel.rxnScores=rxnScores; 0070 tRefModel=simplifyModel(tRefModel,false,false,true); 0071 0072 %This is a rather ugly solution to the issue that it's a bit tricky to keep 0073 %track of which scores belong to which reactions. This requires that 0074 %removeReactions and mergeModels are modified to check for the new field. 0075 tModel.rxnScores=zeros(numel(tModel.rxns),1); 0076 0077 %First merge all models into one big one 0078 fullModel = tRefModel;%mergeModels([{tModel};models],'metNames',true); 0079 0080 %Add that net production is ok 0081 if allowNetProduction==true 0082 %A second column in tModel.b means that the b field is lower and upper 0083 %bound on the RHS. 0084 tModel.b=[tModel.b(:,1) inf(numel(tModel.mets),1)]; 0085 fullModel.b=[fullModel.b(:,1) inf(numel(fullModel.mets),1)]; 0086 end 0087 0088 0089 fullModel.c(:)=0; 0090 0091 %Then minimize for the number of fluxes used. The fixed rxns doesn't need 0092 %to participate 0093 templateRxns = find(~ismember(fullModel.rxns, tModel.rxns)); %Check if this is slow, in that case keep track of this in fitTasksOpt and send it in 0094 0095 [~, J, exitFlag]=ftINITFillGapsMILP(fullModel,templateRxns,params,fullModel.rxnScores(templateRxns),verbose);%only the scores from the template rxns are used, so the others doesn't matter 0096 0097 %Remove everything except for the added ones 0098 addedRxns = fullModel.rxns(templateRxns(J)); 0099 rxnsToAdd.rxns = addedRxns; 0100 rxnsToAdd.equations = constructEquations(fullModel, addedRxns); 0101 rxnsToAdd.ub = fullModel.ub(templateRxns(J)); 0102 rxnsToAdd.lb = fullModel.lb(templateRxns(J)); 0103 newModel = addRxns(origModel,rxnsToAdd, 3, [], true); %we ignore gene rules etc here, they are not needed - we regenerate the model in the end by subtracting rxns from the full model. 0104 0105 end