MATLAB Script
From 2011.igem.org
clear all; close all; clc; changeCobraSolver('glpk','LP'); model = readCbModel('Ec_iJR904_GlcMM.xml');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added reactions %%%%%%%%%%%%%%%%%%%%%%%%%%%%% model = addReaction(model,'Acyl-ACP_Reductase',{'ACP[c]','nadph[c]',...
'Aldehyde','nadp[c]'},[-1,-1,1,1],false,-1000,1000);
model = addReaction(model,'Aldehyde_Decarbonylase',{'Aldehyde',...
'nadph[c]', 'Alkane','nadp[c]','for[c]'},[-1,-1,1,1,1],false,0,1000);
model = addReaction(model,'Alkane_Exchange',{'Alkane','Alkane[e]'},...
[-1,1],false,0,1000);
%%%%%%%%%%%%%%%%%%%%%%%%%% Glucose exchange rate %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% http://www.pnas.org/content/101/8/2235/T3.expansion.html %%%%%%%%% model = changeRxnBounds(model, 'EX_glc(e)', -11.9, 'l');
%%%%%%%%%%%%%%%%%%%%%%%%% Thiamine exchange rate %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% Assumption: exchange rate is proportional to that of glucose %%%%%% model = changeRxnBounds(model, 'THMabc', -0.0004,'l');
% model = changeRxnBounds(model, 'GLYCt', -4, 'l'); % model = changeRxnBounds(model, 'SULabc', -2,'l');
% model = changeObjective(model,{'ACPS1','Alkane_Exchange'}); model = changeObjective(model,{'BiomassEcoli'});
%%%%%%%%%%%%%%%%%%%%%% Removing all existing bounds %%%%%%%%%%%%%%%%%%%%%%% % model = changeRxnBounds(model, 'EX_glc(e)', -1000, 'l'); % model = changeRxnBounds(model, 'EX_o2(e)', -1000, 'l'); % model = changeRxnBounds(model, 'ATPM', -1000, 'l'); % model = changeRxnBounds(model, 'ATPM', 1000, 'u');
%%%%%%%%%%%%%%% Making the synthesis rate of A-ACP nonzero %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% Manual limiting of other A-ACP reactions %%%%%%%%%%%%%%%%% % model = changeRxnBounds(model, 'ACPS1', 10, 'l'); % model = changeRxnBounds(model, 'AACPS1', 10, 'u'); % model = changeRxnBounds(model, 'AACPS2', 10, 'u'); % model = changeRxnBounds(model, 'AACPS3', 10, 'u'); % model = changeRxnBounds(model, 'AACPS4', 10, 'u'); % model = changeRxnBounds(model, 'AACPS5', 10, 'u');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% detectDeadEnds.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%% changeCobraSolver('glpk','MILP'); temp = detectDeadEnds(model); deadEndMets = model.mets(temp); % gapFind(model)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% optGene.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% =%%%%%%%%%%% Requires fuctioning gurobi linear programming solver %%%%%%%%%% % [x,poluation,scores,optGeneSol]=optGene(model,... % 'Aldehyde_Decarbonylase','BiomassEcoli', model.rxns(1:end-3),5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GDLS.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs]=GDLS(model,...
{'ACPS1','BiomassEcoli'});
model2=removeRxns(model,gdlsSolution.KOs);
%%%%%%%%%%%%%%%%%%% Optimize model with biomass function %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Print solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%% solution = optimizeCbModel(model,'max'); % Original model printFluxVector(model, solution.x,true); solution2 = optimizeCbModel(model2,'max'); % Model with reaction knock outs printFluxVector(model2, solution2.x,true);