% This is the solution of excersice 2

%% Setting the Matlab enviroment
clear;
clc;
restoredefaultpath;

addpath(genpath('/Users/davidliaskos/GIT_Folders/mattfa'))
addpath(genpath('/Users/davidliaskos/IBM/CPLEX_Studio1210/cplex/matlab'))
changeCobraSolver('cplex_direct', 'LP')
addpath(genpath('./'))

% Load E. coli model
load('small_ecoli.mat');

%% Part 1: Simulating Growth

% FBA aerobic
excRxns={'DM_glc_e','DM_lac-D_e','DM_ac_e','DM_etoh_e'};
Substrates={'D-glucose','D-lactate','Acetate','Ethanol'};
model=model_red;

model.lb(find(ismember(model.rxns, 'DM_o2_e'))) = -20;
model.ub(find(ismember(model.rxns, 'DM_o2_e'))) = 0;

maxYield_table{1,1}='Substrate';
maxYield_table{1,2}='Max. yield aerobic FBA';
maxYield_table{1,3}='Max. yield aerobic TFA ';
maxYield_table{1,4}='Max. yield anaerobic FBA';
maxYield_table{1,5}='Max. yield anaerobic TFA ';
 
for i=2:size(excRxns,2)+1
    maxYield_table{i,1}=Substrates{i-1};
end
 
for i=1:length(excRxns)
     
    modeli=model;
    
    carbonSource=excRxns{i};
    excRxns_rest=excRxns(find(not(ismember(excRxns,carbonSource))));

    modeli.lb(find(ismember(modeli.rxns, carbonSource))) = -10;
    modeli.lb(find(ismember(modeli.rxns, excRxns_rest))) = 0;
    
    FBAsolutioni=solveFBAmodelCplex(modeli);
    maxYield_table{i+1,2}=FBAsolutioni.f;
end

% FBA anaerobic
model.lb(find(ismember(model.rxns, 'DM_o2_e'))) = 0;

for i=1:length(excRxns)
     
    modeli=model;
    
    carbonSource=excRxns{i};
    excRxns_rest=excRxns(find(not(ismember(excRxns,carbonSource))));

    modeli.lb(find(ismember(modeli.rxns, carbonSource))) = -10;
    modeli.lb(find(ismember(modeli.rxns, excRxns_rest))) = 0;
    
    FBAsolutioni=solveFBAmodelCplex(modeli);
    maxYield_table{i+1,4}=FBAsolutioni.f;
end

% Convert model to TFA
load('thermo_data.mat')

model_tmp = prepModelforTFA(model_red, DB_AlbertyUpdate,model_red.CompartmentData);
[tmodel] = convToTFA(model_tmp, DB_AlbertyUpdate, {}, 'DGo', [], 0.0);
tmodel = addNetFluxVariables(tmodel);

% TFA aerobic
model=tmodel;

model.var_ub(find(ismember(model.varNames, 'R_DM_o2_e'))) = 20;

for i=1:length(excRxns)
     
    modeli=model;
    
    carbonSource=excRxns{i};
    excRxns_rest=excRxns(find(not(ismember(excRxns,carbonSource))));

    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',carbonSource)))) = 10;
    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',excRxns_rest)))) = 0;
    
    TFAsolutioni=solveTFAmodelCplex(modeli);
    maxYield_table{i+1,3}=TFAsolutioni.val;
end

% TFA anaerobic
model=tmodel;

model.var_ub(find(ismember(model.varNames, 'R_DM_o2_e'))) = 0;

for i=1:length(excRxns)
     
    modeli=model;
    
    carbonSource=excRxns{i};
    excRxns_rest=excRxns(find(not(ismember(excRxns,carbonSource))));

    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',carbonSource)))) = 10;
    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',excRxns_rest)))) = 0;
    
    TFAsolutioni=solveTFAmodelCplex(modeli);
    maxYield_table{i+1,5}=TFAsolutioni.val;
end

%% Part 2: Integrating concentration data
metabolomics = readtable('metabolomics_data.csv');

metabolomics(:,1)=strcat('LC_',metabolomics{:,1});

for i=1:size(metabolomics,1)
    if ~isempty(find(ismember(tmodel.varNames,metabolomics{i,1}))) 
        if ~isnan(metabolomics{i,4})
            
            if (metabolomics{i,3}-metabolomics{i,4}) > 0
                tmodel.var_ub(find(ismember(tmodel.varNames,metabolomics{i,1})))=log(1e-3*(metabolomics{i,3}+metabolomics{i,4}));
                tmodel.var_lb(find(ismember(tmodel.varNames,metabolomics{i,1})))=log(1e-3*(metabolomics{i,3}-metabolomics{i,4}));
            elseif ((metabolomics{i,3}-metabolomics{i,4}) <= 0) && (tmodel.var_lb(find(ismember(tmodel.varNames,metabolomics{i,1}))) <= log(1e-3*(metabolomics{i,3}+metabolomics{i,4})))
                tmodel.var_ub(find(ismember(tmodel.varNames,metabolomics{i,1})))=log(1e-3*(metabolomics{i,3}+metabolomics{i,4}));
            end

        else

            tmodel.var_ub(find(ismember(tmodel.varNames,metabolomics{i,1})))=log(1e-3*(metabolomics{i,3}*1.5));
            tmodel.var_lb(find(ismember(tmodel.varNames,metabolomics{i,1})))=log(1e-3*(metabolomics{i,3}*0.5));

        end
    end
end

 maxYield_table_c{1,1}='Substrate';
 maxYield_table_c{1,2}='Max. yield aerobic TFA metabolomic';
 maxYield_table_c{1,3}='Max. yield anaerobic TFA metabolimics ';
 for i=2:size(excRxns,2)+1
    maxYield_table_c{i,1}=Substrates{i-1};
 end

 % TFA aerobic
model=tmodel;

model.var_ub(find(ismember(model.varNames, 'R_DM_o2_e'))) = 20;

for i=1:length(excRxns)
     
    modeli=model;
    
    carbonSource=excRxns{i};
    excRxns_rest=excRxns(find(not(ismember(excRxns,carbonSource))));

    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',carbonSource)))) = 10;
    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',excRxns_rest)))) = 0;
    
    TFAsolutioni=solveTFAmodelCplex(modeli);
    maxYield_table_c{i+1,2}=TFAsolutioni.val;
end

% TFA anaerobic
model=tmodel;

model.var_ub(find(ismember(model.varNames, 'R_DM_o2_e'))) = 0;

for i=1:length(excRxns)
     
    modeli=model;
    
    carbonSource=excRxns{i};
    excRxns_rest=excRxns(find(not(ismember(excRxns,carbonSource))));

    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',carbonSource)))) = 10;
    modeli.var_ub(find(ismember(modeli.varNames, strcat('R_',excRxns_rest)))) = 0;
    
    TFAsolutioni=solveTFAmodelCplex(modeli);
    maxYield_table_c{i+1,3}=TFAsolutioni.val;
end

save('./WS_exc2')