Team:IISER-Pune-India/TeamNotebook/Team Notebook 11e9906b405543a29361e38a50c81024/E coli Model 35df109f96a148ea8c27d813ca507dfe

E.coli Model

E.coli Model

Date
DepartmentDry Lab
DescriptionBottom-most code is the final model
HP sub-branch
Links/media
Participants
Property
Property 1
Property 2
TypeDeliverables

Here model is the iML1515 model

%Version 1 (update available)
filename = 'iML1515.mat'
model = readCbModel(filename)

%Knockouts

%pta
model = changeRxnBounds(model,'PTAr',0,'b')
model = changeRxnBounds(model,'PTA2',0,'b')

%adhE
model = changeRxnBounds(model,'ACALD',0,'b')
model = changeRxnBounds(model,'ALCD2x',0,'b')
model = changeRxnBounds(model,'ALCD19',0,'b')

%frdA
model = changeRxnBounds(model,'FRD2',0,'b')
model = changeRxnBounds(model,'FRD3',0,'b')

%ldhA
model = changeRxnBounds(model,'LDH_D',0,'b')


%Heterologous Reactions

%Clostridium acetobutylicum

%(S)-3-hydroxybutanoyl-CoA + NAD+ ↔ acetoacetyl-CoA + NADH + H+ 
model = addReaction(model, 'HBD','metaboliteList', {'3hbcoa_c', 'nad_c', 'aacoa_c', 'nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'hbd'}, 'checkDuplicate', true)  
disp('HBD done')

%(S)-3-hydroxyhexanoyl-CoA → (2E)-hexenoyl-CoA + H2O [already exists]
model = addReaction(model, 'CRT_1','metaboliteList', {'3hhcoa_c', 'hx2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', false,'geneNameList', {'crt'}, 'checkDuplicate', true) 
disp('CRT_1 done') %ECOAH2 already exists, should we overexpress?

%(S)-3-hydroxybutanoyl-CoA ↔ crotonyl-CoA + H2O
model = addReaction(model, 'CRT_2','metaboliteList', {'3hhcoa_c', 'b2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible',false, 'geneNameList', {'crt'}, 'checkDuplicate', true)
disp('CRT_2 done')

%1-hexanal + NADH + H+ → hexan-1-ol + NAD+
model = addReaction(model, 'ADHE2_1','metaboliteList', {'1_hexal_c', 'nadh_c', 'h_c', '1_hexol_c', 'nad_c' },'stoichCoeffList', [-1; -1; -1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)  % 1_hexal_c - 1-hexanal, 1_hexol_c - hexan-1-ol
disp('ADHE2_1 done')

%hexanoyl-CoA + NAD(P)H + H+ → 1-hexanal + coenzyme A + NAD(P)+
model = addReaction(model, 'ADHE2_2','metaboliteList', {'hxcoa_c', 'nadph_c', 'h_c', '1_hexal_c', 'coa_c', 'nadp_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_2 done')

%butan-1-ol + NAD+ ↔ 1-butanal + NADH + H+
model = addReaction(model, 'ADHE2_3','metaboliteList', {'1btol_e','nad_c', 'btal_c','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true) % 1btol_e - butanol
disp('ADHE2_3 done')

%1-butanal + coenzyme A + NAD(P)+ ↔ butanoyl-CoA + NAD(P)H + H+
model = addReaction(model, 'ADHE2_4','metaboliteList', {'btal_c','coa_c', 'nadp_c','btcoa_c','nadph_c', 'h_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_4 done')

%Treponema denticola
%acyl-CoA + NAD+ = trans-2,3-dehydroacyl-CoA + NADH + H+
model = addReaction(model, 'TER','metaboliteList', {'acoa_c' 'nad_c','t2ecoa','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'ter'}, 'checkDuplicate', true)  %t2ecoa - trans-2,3-dehydroacyl-CoA 
disp('TER done')

%add reactions from csc genes
model = addReaction(model, 'CSC_B', 'metaboliteList', {'sucr_p', 'h_p', 'sucr_c', 'h_c'}, 'stoichCoeffList', [-1;-1;1;1]);
model = addReaction(model, 'CSC_A', 'metaboliteList', {'sucr_c', 'h2o_c', 'glc__D_c', 'fru_c'}, 'stoichCoeffList', [-1;-1;1;1], 'reversible', false);
model = addReaction(model, 'CSC_K', 'metaboliteList', {'atp_c', 'fru_c', 'h_c', 'adp_c', 'f6p_c'}, 'stoichCoeffList', [1;1;-1;-1;-1], 'reversible', false);

mets_length_f = length(model.mets)
rxns_length_f = length(model.rxns)

%Adding genes to the appropriate reactions

model = addGenes(model, {'hbd', 'crt', 'adhE2', 'ter'})

Rxns = model.rxns(end-6:end)
Genes = model.genes(end-3:end)

model.rxnGeneMat(end-6, end-3)= 1
model.rxnGeneMat(end-5, end-2)= 1
model.rxnGeneMat(end-4:end-1 , end-1)= 1
model.rxnGeneMat(end , end)= 1

%First optimization for biomass

changeCobraSolver('gurobi', 'all')
solution_1 = optimizeCbModel(model)
%Version 2 (update available)
initCobraToolbox

model = readCbModel('ecoli.mat')

mets_length_i = length(model.mets)
rxns_length_i = length(model.rxns)

%pta
model = changeRxnBounds(model,'PTAr',0,'b')
model = changeRxnBounds(model,'PTA2',0,'b')

%adhE
model = changeRxnBounds(model,'ACALD',0,'b')
model = changeRxnBounds(model,'ALCD2x',0,'b')
model = changeRxnBounds(model,'ALCD19',0,'b')

%frdA
model = changeRxnBounds(model,'FRD2',0,'b')
model = changeRxnBounds(model,'FRD3',0,'b')

%ldhA
model = changeRxnBounds(model,'LDH_D',0,'b')

%anaerobic
model = changeRxnBounds(model,'EX_o2_e',0,'b')

%Clostridium acetobutylicum

%(S)-3-hydroxybutanoyl-CoA + NAD+ ↔ acetoacetyl-CoA + NADH + H+ 
model = addReaction(model, 'HBD','metaboliteList', {'3hbcoa_c', 'nad_c', 'aacoa_c', 'nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'hbd'}, 'checkDuplicate', true)  
disp('HBD done')

%(S)-3-hydroxyhexanoyl-CoA → (2E)-hexenoyl-CoA + H2O [already exists]
model = addReaction(model, 'CRT_1','metaboliteList', {'3hhcoa_c', 'hx2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', false,'geneNameList', {'crt'}, 'checkDuplicate', true) 
disp('CRT_1 done') %ECOAH2 already exists, should we overexpress?

%(S)-3-hydroxybutanoyl-CoA ↔ crotonyl-CoA + H2O  (change reversibility)
model = addReaction(model, 'CRT_2','metaboliteList', {'3hhcoa_c', 'b2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', true, 'geneNameList', {'crt'}, 'checkDuplicate', true)
disp('CRT_2 done') 

%1-hexanal + NADH + H+ → hexan-1-ol + NAD+
model = addReaction(model, 'ADHE2_1','metaboliteList', {'1_hexal_c', 'nadh_c', 'h_c', '1_hexol_c', 'nad_c' },'stoichCoeffList', [-1; -1; -1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)  % 1_hexal_c - 1-hexanal, 1_hexol_c - hexan-1-ol
disp('ADHE2_1 done')

%hexanoyl-CoA + NAD(P)H + H+ → 1-hexanal + coenzyme A + NAD(P)+
model = addReaction(model, 'ADHE2_2','metaboliteList', {'hxcoa_c', 'nadph_c', 'h_c', '1_hexal_c', 'coa_c', 'nadp_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_2 done')

%butan-1-ol + NAD+ ↔ 1-butanal + NADH + H+
model = addReaction(model, 'ADHE2_3','metaboliteList', {'1btol_c','nad_c', 'btal_c','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true) % 1btol_e - butanol
disp('ADHE2_3 done')

%1-butanal + coenzyme A + NAD(P)+ ↔ butanoyl-CoA + NAD(P)H + H+
model = addReaction(model, 'ADHE2_4','metaboliteList', {'btal_c','coa_c', 'nadp_c','btcoa_c','nadph_c', 'h_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_4 done')

%Treponema denticola
%acyl-CoA + NAD+ = trans-2,3-dehydroacyl-CoA + NADH + H+
model = addReaction(model, 'TER','metaboliteList', {'acoa_c' 'nad_c','t2ecoa','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'ter'}, 'checkDuplicate', true)  %t2ecoa - trans-2,3-dehydroacyl-CoA 
disp('TER done')

mets_length_f = length(model.mets)
rxns_length_f = length(model.rxns)
%disp('initial_mets',mets_length_i, 'final_mets', mets_length_f, 'initial_rxns', rxns_length_i,'final_rxns', rxns_length_f)

model = addReaction(model,'EX_1btol_e[e]','metaboliteList', {'1btol_e'} ,'stoichCoeffList', [-1])    %butanol exchange reaction

model = addGenes(model, {'hbd', 'crt', 'adhE2', 'ter'})

Rxns = model.rxns(end-7:end)
Genes = model.genes(end-3:end)

model.rxnGeneMat(end-7, end-3)= 1
model.rxnGeneMat(end-6, end-2)= 1
model.rxnGeneMat(end-5:end-2 , end-1)= 1
model.rxnGeneMat(end-1, end)= 1

changeCobraSolver('gurobi', 'all')
solution_1 = optimizeCbModel(model)

solution_1.v(end-6:end)
solution_1.v(255)           %CRT_1 reaction to be overexpressed
solution_1.v(1537)          %atoB reaction to be overexpressed5
%Version 3
initCobraToolbox

model = readCbModel('ecoli.mat')

mets_length_i = length(model.mets)
rxns_length_i = length(model.rxns)

%pta
model = changeRxnBounds(model,'PTAr',0,'b')
model = changeRxnBounds(model,'PTA2',0,'b')

%adhE
model = changeRxnBounds(model,'ACALD',0,'b')
model = changeRxnBounds(model,'ALCD2x',0,'b')
model = changeRxnBounds(model,'ALCD19',0,'b')

%frdA
model = changeRxnBounds(model,'FRD2',0,'b')
model = changeRxnBounds(model,'FRD3',0,'b')

%ldhA
model = changeRxnBounds(model,'LDH_D',0,'b')

%anaerobic
model = changeRxnBounds(model,'EX_o2_e',0,'b')

%sucrose exchange
model.lb(298) = -10 

%no glucose uptake, 'EX_glc__D_e' ID - 181
model = changeRxnBounds(model,'EX_glc__D_e',0,'b')

%Clostridium acetobutylicum

%(S)-3-hydroxybutanoyl-CoA + NAD+ ↔ acetoacetyl-CoA + NADH + H+ 
model = addReaction(model, 'HBD','metaboliteList', {'3hbcoa_c', 'nad_c', 'aacoa_c', 'nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'hbd'}, 'checkDuplicate', true)  
disp('HBD done')

%(S)-3-hydroxyhexanoyl-CoA → (2E)-hexenoyl-CoA + H2O [already exists]
model = addReaction(model, 'CRT_1','metaboliteList', {'3hhcoa_c', 'hx2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', false,'geneNameList', {'crt'}, 'checkDuplicate', true) 
disp('CRT_1 done') %ECOAH2 already exists, should we overexpress?

%(S)-3-hydroxybutanoyl-CoA → crotonyl-CoA + H2O 
model = addReaction(model, 'CRT_2','metaboliteList', {'3hhcoa_c', 'b2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', false, 'geneNameList', {'crt'}, 'checkDuplicate', true)
disp('CRT_2 done') 

%1-hexanal + NADH + H+ → hexan-1-ol + NAD+
model = addReaction(model, 'ADHE2_1','metaboliteList', {'1_hexal_c', 'nadh_c', 'h_c', '1_hexol_c', 'nad_c' },'stoichCoeffList', [-1; -1; -1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)  % 1_hexal_c - 1-hexanal, 1_hexol_c - hexan-1-ol
disp('ADHE2_1 done')

%hexanoyl-CoA + NAD(P)H + H+ → 1-hexanal + coenzyme A + NAD(P)+
model = addReaction(model, 'ADHE2_2','metaboliteList', {'hxcoa_c', 'nadph_c', 'h_c', '1_hexal_c', 'coa_c', 'nadp_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_2 done')

%butan-1-ol + NAD+ ↔ 1-butanal + NADH + H+
model = addReaction(model, 'ADHE2_3','metaboliteList', {'1btol_c','nad_c', 'btal_c','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true) % 1btol_e - butanol
disp('ADHE2_3 done')

%1-butanal + coenzyme A + NAD(P)+ ↔ butanoyl-CoA + NAD(P)H + H+
model = addReaction(model, 'ADHE2_4','metaboliteList', {'btal_c','coa_c', 'nadp_c','btcoa_c','nadph_c', 'h_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_4 done')

%butan-1-ol ↔
model = addReaction(model, '1btol_trnpt', 'metaboliteList', {'1btol_c', '1btol_e'}, 'stoichCoeffList', [-1; 1])
disp('1btol_trnpt done')

%Treponema denticola
%acyl-CoA + NAD+ = trans-2,3-dehydroacyl-CoA + NADH + H+
model = addReaction(model, 'TER','metaboliteList', {'acoa_c' 'nad_c','t2ecoa','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'ter'}, 'checkDuplicate', true)  %t2ecoa - trans-2,3-dehydroacyl-CoA 
disp('TER done')

%csc reactions
model = addReaction(model, 'CSC_B', 'metaboliteList', {'sucr_p', 'h_p', 'sucr_c', 'h_c'}, 'stoichCoeffList', [-1;-1;1;1])
model = addReaction(model, 'CSC_A', 'metaboliteList', {'sucr_c', 'h2o_c', 'glc__D_c', 'fru_c'}, 'stoichCoeffList', [-1;-1;1;1], 'reversible', false)
model = addReaction(model, 'CSC_K', 'metaboliteList', {'atp_c', 'fru_c', 'h_c', 'adp_c', 'f6p_c'}, 'stoichCoeffList', [1;1;-1;-1;-1], 'reversible', false)

mets_length_f = length(model.mets)
rxns_length_f = length(model.rxns)
%disp('initial_mets',mets_length_i, 'final_mets', mets_length_f, 'initial_rxns', rxns_length_i,'final_rxns', rxns_length_f)

model = addReaction(model,'EX_1btol_e','metaboliteList', {'1btol_e'} ,'stoichCoeffList', [-1])    %butanol exchange reaction

model = addGenes(model, {'hbd', 'crt', 'adhE2', 'ter', 'csc'})

Rxns = model.rxns(end-11:end)
Genes = model.genes(end-4:end)

model.rxnGeneMat(2713, 1517)= 1
model.rxnGeneMat(2714, 1518)= 1
model.rxnGeneMat(2715:2718, 1519)= 1
model.rxnGeneMat(2720, 1520)= 1
model.rxnGeneMat(2721:2723, 1521)= 1

model = changeRxnBounds(model,'EX_1btol_e',0,'l')
model = changeRxnBounds(model,'1btol_trnpt',0,'l')
%Version 4
initCobraToolbox

model = readCbModel('ecoli.mat')

mets_length_i = length(model.mets)
rxns_length_i = length(model.rxns)

%pta
model = changeRxnBounds(model,'PTAr',0,'b')
model = changeRxnBounds(model,'PTA2',0,'b')

%adhE
model = changeRxnBounds(model,'ACALD',0,'b')
model = changeRxnBounds(model,'ALCD2x',0,'b')
model = changeRxnBounds(model,'ALCD19',0,'b')

%frdA
model = changeRxnBounds(model,'FRD2',0,'b')
model = changeRxnBounds(model,'FRD3',0,'b')

%ldhA
model = changeRxnBounds(model,'LDH_D',0,'b')

%anaerobic
model = changeRxnBounds(model,'EX_o2_e',0,'b')

%sucrose exchange
model.lb(298) = -10 

%no glucose uptake, 'EX_glc__D_e' ID - 181
model = changeRxnBounds(model,'EX_glc__D_e',0,'b')

%Clostridium acetobutylicum

%(S)-3-hydroxybutanoyl-CoA + NAD+ ↔ acetoacetyl-CoA + NADH + H+ 
model = addReaction(model, 'HBD','metaboliteList', {'3hbcoa_c', 'nad_c', 'aacoa_c', 'nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'hbd'}, 'checkDuplicate', true)  
disp('HBD done')

%(S)-3-hydroxyhexanoyl-CoA → (2E)-hexenoyl-CoA + H2O [already exists: 'ECOAH2', 255]
model = addReaction(model, 'CRT_1','metaboliteList', {'3hhcoa_c', 'hx2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', false,'geneNameList', {'crt'}, 'checkDuplicate', true) 
disp('CRT_1 done') %ECOAH2 already exists, should we overexpress?

%(S)-3-hydroxybutanoyl-CoA ↔ crotonyl-CoA + H2O [already exists: 'ECOAH1', 254] (reversibility changed)
model = addReaction(model, 'CRT_2','metaboliteList', {'3hbcoa_c', 'b2coa_c', 'h2o_c'},'stoichCoeffList', [-1; 1; 1], 'reversible', false, 'geneNameList', {'crt'}, 'checkDuplicate', true)
disp('CRT_2 done') 

%1-hexanal + NADH + H+ → hexan-1-ol + NAD+
model = addReaction(model, 'ADHE2_1','metaboliteList', {'1_hexal_c', 'nadh_c', 'h_c', '1_hexol_c', 'nad_c' },'stoichCoeffList', [-1; -1; -1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)  % 1_hexal_c - 1-hexanal, 1_hexol_c - hexan-1-ol
disp('ADHE2_1 done')

%hexanoyl-CoA + NAD(P)H + H+ → 1-hexanal + coenzyme A + NAD(P)+
model = addReaction(model, 'ADHE2_2','metaboliteList', {'hxcoa_c', 'nadph_c', 'h_c', '1_hexal_c', 'coa_c', 'nadp_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_2 done')

%butan-1-ol + NAD+ ↔ 1-butanal + NADH + H+
model = addReaction(model, 'ADHE2_3','metaboliteList', {'1btol_c','nad_c', 'btal_c','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true) % 1btol_e - butanol
disp('ADHE2_3 done')

%1-butanal + coenzyme A + NAD(P)+ ↔ butanoyl-CoA + NAD(P)H + H+
model = addReaction(model, 'ADHE2_4','metaboliteList', {'btal_c','coa_c', 'nadp_c','btcoa_c','nadph_c', 'h_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', true)
disp('ADHE2_4 done')

%butan-1-ol_c → butanol-1-ol_e
model = addReaction(model, '1btol_trnpt', 'metaboliteList', {'1btol_c', '1btol_e'}, 'stoichCoeffList', [-1; 1])
disp('1btol_trnpt done')

%Treponema denticola
%acyl-CoA + NAD+ = trans-2,3-dehydroacyl-CoA + NADH + H+
model = addReaction(model, 'TER','metaboliteList', {'acoa_c' 'nad_c','t2ecoa','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'ter'}, 'checkDuplicate', true)  %t2ecoa - trans-2,3-dehydroacyl-CoA 
disp('TER done')

%csc reactions
model = addReaction(model, 'CSC_B', 'metaboliteList', {'sucr_p', 'h_p', 'sucr_c', 'h_c'}, 'stoichCoeffList', [-1;-1;1;1], 'checkDuplicate', true)
model = addReaction(model, 'CSC_A', 'metaboliteList', {'sucr_c', 'h2o_c', 'glc__D_c', 'fru_c'}, 'stoichCoeffList', [-1;-1;1;1], 'reversible', false, 'checkDuplicate', true)
model = addReaction(model, 'CSC_K', 'metaboliteList', {'atp_c', 'fru_c', 'h_c', 'adp_c', 'f6p_c'}, 'stoichCoeffList', [1;1;-1;-1;-1], 'reversible', false, 'checkDuplicate', true)

mets_length_f = length(model.mets)
rxns_length_f = length(model.rxns)
%disp('initial_mets',mets_length_i, 'final_mets', mets_length_f, 'initial_rxns', rxns_length_i,'final_rxns', rxns_length_f)

model = addReaction(model,'EX_1btol_e','metaboliteList', {'1btol_e'} ,'stoichCoeffList', [-1])    %butanol exchange reaction

model = addGenes(model, {'hbd', 'adhE2', 'ter', 'csc'})

Rxns = model.rxns(end-10:end)
Genes = model.genes(end-3:end)

model.rxnGeneMat(2713, 1517)= 1
model.rxnGeneMat(2714:2717, 1518)= 1
model.rxnGeneMat(2719, 1519)= 1
model.rxnGeneMat(2720:2722, 1520)= 1

model = changeRxnBounds(model,'EX_1btol_e',0,'l')
model = changeRxnBounds(model,'1btol_trnpt',0,'l')


changeCobraSolver('gurobi', 'all')
solution_1 = optimizeCbModel(model)

solution_1.v(2713:2723)
solution_1.v(254)           %CRT_2 reaction to be overexpressed
solution_1.v(255)           %CRT_1 reaction to be overexpressed
solution_1.v(1537)          %atoB reaction to be overexpressed

model = changeRxnBounds(model,'CSC_K',10,'u')  %testing
changeCobraSolver('gurobi', 'all')
solution_1 = optimizeCbModel(model)

solution_1.v(2713:2723)
solution_1.v(254)           %CRT_2 reaction to be overexpressed
solution_1.v(255)           %CRT_1 reaction to be overexpressed
solution_1.v(1537)          %atoB reaction to be overexpressed

changeCobraSolver('gurobi', 'all')
counter=0
flux_CRT_1 = nan(125,1)
flux_CRT_2 = nan(125,1)
flux_atob = nan(125,1)
butanol_exc = nan(125,1)
growth = nan(125,1)
for i = 12:0.2:13
    model = changeRxnBounds(model,'ECOAH1',i,'b')
    for j = 12:0.2:13
        model = changeRxnBounds(model,'ECOAH2',j,'b')
        for k = 12:0.2:13
            model= changeRxnBounds(model, 'ACACT1r', k, 'b')
            opt = optimizeCbModel(model)
            counter=counter+1
            if opt.stat==1
                flux_CRT_2(counter)= i
                flux_CRT_1(counter)= j
                flux_atob(counter) = k
                butanol_exc(counter)= opt.v(2723)
                growth(counter) = opt.f
            else
                flux_CRT_1(counter)= i
                flux_CRT_2(counter)=j
                flux_atob(counter) = k
                butanol_exc(counter)= 0
                growth(counter) = 0
            end    
        end
    end
end

grph = table(flux_CRT_1, flux_CRT_2, flux_atob, butanol_exc, growth)
writetable(grph,'ecoli_ovrxps_fluxes.xlsx')

scatter3(grph.flux_CRT_1, grph.flux_CRT_2,  grph.flux_atob,'S', grph.butanol_exc,'C', grph.growth, 'filled')
ax = gca;
xlabel('Flux through crt_1 gene product')
ylabel('Flux through crt_2 gene product')
zlabel('Flux through atoB gene product')

cb = colorbar;                                     % create and label the colorbar
cb.Label.String = 'Growth rate';
savefig('Ecoli_butnl_grwth.fig')
saveas(gcf, 'Ecoli_butnl_grwth.png')
%Final E coli model
initCobraToolbox

model = readCbModel('ecoli.mat')

%pta
model = changeRxnBounds(model,'PTAr',0,'b')
model = changeRxnBounds(model,'PTA2',0,'b')

%adhE
model = changeRxnBounds(model,'ACALD',0,'b')
model = changeRxnBounds(model,'ALCD2x',0,'b')
model = changeRxnBounds(model,'ALCD19',0,'b')

%frdA
model = changeRxnBounds(model,'FRD2',0,'b')
model = changeRxnBounds(model,'FRD3',0,'b')

%ldhA
model = changeRxnBounds(model,'LDH_D',0,'b')

%anaerobic
model = changeRxnBounds(model,'EX_o2_e',0,'b')

%sucrose exchange
model.lb(298) = -10 

%no glucose uptake, 'EX_glc__D_e' ID - 181
model = changeRxnBounds(model,'EX_glc__D_e',0,'b')

%Clostridium acetobutylicum 

%1-hexanal + NADH + H+ → hexan-1-ol + NAD+
model = addReaction(model, 'ADHE2_1','metaboliteList', {'1_hexal_c', 'nadh_c', 'h_c', '1_hexol_c', 'nad_c' },'stoichCoeffList', [-1; -1; -1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', 1);  % 1_hexal_c - 1-hexanal, 1_hexol_c - hexan-1-ol
disp('ADHE2_1 done')

%hexanoyl-CoA + NAD(P)H + H+ → 1-hexanal + coenzyme A + NAD(P)+
model = addReaction(model, 'ADHE2_2','metaboliteList', {'hxcoa_c', 'nadph_c', 'h_c', '1_hexal_c', 'coa_c', 'nadp_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', false, 'geneNameList', {'adhE2'}, 'checkDuplicate', 1);
disp('ADHE2_2 done')

%butan-1-ol + NAD+ ↔ 1-butanal + NADH + H+
model = addReaction(model, 'ADHE2_3','metaboliteList', {'1btol_c','nad_c', 'btal_c','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', 1); % 1btol_e - butanol
disp('ADHE2_3 done')

%1-butanal + coenzyme A + NAD(P)+ ↔ butanoyl-CoA + NAD(P)H + H+
model = addReaction(model, 'ADHE2_4','metaboliteList', {'btal_c','coa_c', 'nadp_c','btcoa_c','nadph_c', 'h_c' },'stoichCoeffList', [-1; -1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'adhE2'}, 'checkDuplicate', 1);
disp('ADHE2_4 done')

%butan-1-ol_c → butanol-1-ol_e
model = addReaction(model, '1btol_trnpt', 'metaboliteList', {'1btol_c', '1btol_e'}, 'stoichCoeffList', [-1; 1]);
disp('1btol_trnpt done')

%Treponema denticola
%acyl-CoA + NAD+ = trans-2,3-dehydroacyl-CoA + NADH + H+
model = addReaction(model, 'TER','metaboliteList', {'acoa_c' 'nad_c','t2ecoa','nadh_c', 'h_c' },'stoichCoeffList', [-1; -1; 1; 1; 1], 'reversible', true, 'geneNameList', {'ter'}, 'checkDuplicate', 1);  %t2ecoa - trans-2,3-dehydroacyl-CoA 
disp('TER done')

%csc reactions
model = addReaction(model, 'CSC_B', 'metaboliteList', {'sucr_p', 'h_p', 'sucr_c', 'h_c'}, 'stoichCoeffList', [-1;-1;1;1], 'checkDuplicate', 1)
model = addReaction(model, 'CSC_A', 'metaboliteList', {'sucr_c', 'h2o_c', 'glc__D_c', 'fru_c'}, 'stoichCoeffList', [-1;-1;1;1], 'reversible', false, 'checkDuplicate', 1)


model = addReaction(model,'EX_1btol_e','metaboliteList', {'1btol_e'} ,'stoichCoeffList', [-1])    %butanol exchange reaction

model = addGenes(model, {'hbd', 'adhE2', 'ter', 'csc'})

Rxns = model.rxns(2713:2721)
Genes = model.genes(1517: 1520)

model.rxnGeneMat(2713: 2716, 1518)= 1;
model.rxnGeneMat(2718, 1519)= 1;
model.rxnGeneMat(2719:2720, 1520)= 1;

model = changeRxnBounds(model,'EX_1btol_e',0,'l');
model = changeRxnBounds(model,'1btol_trnpt',0,'l');

model = changeRxnBounds(model,{'HACD1', 'ECOAH1', 'ACACT1r'},14,'l');  %hbd, crt_2, atoB
model = changeRxnBounds(model,'ECOAH2',12,'l')  %crt_1
model = changeRxnBounds(model,'HEX7',20,'l')
model = changeRxnBounds(model,'BIOMASS_Ec_iML1515_core_75p37M',0.136,'l')