# 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);