% This is the solution of excersice 1

%% Setting the Matlab enviroment
% Adapt the pathnames based on the locations of your mattfa and CPLEX
% folders location
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('./'))


%% Part 1: Load and compare GEMs
load('gem_ecoli_curated.mat')
load('Yeast.mat')


%% Part 2: Simulate growth

% Find the rxns indices for the uptake of glucose and oxygen
id_glc = find(ismember(tmodel.rxns, 'DM_glc_e')); % Find the index of the glucose exchange reaction
id_o2 = find(ismember(tmodel.rxns, 'DM_o2_e')); % Find the index of the oxygen exchange reaction

% a) Simulate aerobic growth of E. coli

% Define uptake rate of glucose by modyfing the bounds of the reaction
tmodel.lb(id_glc) = -10; % Set the lower bound to -10 mmol/gDW/hr
tmodel.ub(id_glc) = 0; % Dont allow glucose secretion

% Define uptake rate of oxygen by modyfing the bounds of the reaction
tmodel.lb(id_o2) = -100;  % Set the lower bound to -100 mmol/gDW/hr
tmodel.ub(id_o2) = 0; % Dont allow oxygen secretion

% Find the biomass reaction index
id_biomass = find(ismember(tmodel.rxns, 'Ec_biomass_iJO1366_WT_53p95M')); % Find the biomass reaciton index
% Allow Ecoli model not to grow
tmodel.lb(id_biomass) = 0;

% Set the objective function (c matrix) to the biomass reaction
tmodel.c = tmodel.c*0;
tmodel.c(id_biomass) = 1;

% Use the mattfa solveFBAmodelCplex function to optimize the FBA problem
sol = solveFBAmodelCplex(tmodel);   % By default solveFBAmodelCplex will solve a maximization problem
sol.f % Displat the solution
solutions.ecoli_aerobic = sol.f;

% b) Simulate anaerobic growth of E. coli

% Define uptake rate of oxygen by modyfing the bounds of the reaction
tmodel.lb(id_o2) = 0;  % Dont allow oxygen uptake
tmodel.ub(id_o2) = 0; % Dont allow oxygen secretion

% Use the mattfa solveFBAmodelCplex function to optimize the FBA problem
sol = solveFBAmodelCplex(tmodel);
sol.f % Displat the solution
solutions.ecoli_anaerobic = sol.f;


% Perform th same analysis for the yeast model

% Find the rxns indices for the uptake of glucose and oxygen for Yeast
id_glc_yeast = find(ismember(iMM904.rxns, 'EX_glc(e)')); % Find the index of the glucose exchange reaction
id_o2_yeast = find(ismember(iMM904.rxns, 'EX_o2(e)')); % Find the index of the oxygen exchange reaction

% c) Simulate aerobic growth of Yeast

% Define uptake rate of glucose by modyfing the bounds of the reaction
iMM904.lb(id_glc_yeast) = -10; % Set the lower bound to -10 mmol/gDW/hr
iMM904.ub(id_glc_yeast) = 0; % Dont allow glucose secretion

% Define uptake rate of oxygen by modyfing the bounds of the reaction
iMM904.lb(id_o2_yeast) = -100;  % Set the lower bound to -100 mmol/gDW/hr
iMM904.ub(id_o2_yeast) = 0; % Dont allow oxygen secretion

% Find the biomass reaction index
id_biomass_yeast = find(ismember(iMM904.rxns, 'biomass_SC5_notrace')); % Find the biomass reaciton index
% Allow Yeast model not to grow
iMM904.lb(id_biomass_yeast) = 0;

% Set the objective function (c matrix) to the biomass reaction
iMM904.c = iMM904.c*0;
iMM904.c(id_biomass_yeast) = 1;

% Use the mattfa solveFBAmodelCplex function to optimize the FBA problem
sol = solveFBAmodelCplex(iMM904);
sol.f % Displat the solution
solutions.yeast_aerobic = sol.f;

% d) Simulate anaerobic growth of E. coli

% Define uptake rate of oxygen by modyfing the bounds of the reaction
iMM904.lb(id_o2_yeast) = 0;  % Dont allow oxygen uptake
iMM904.ub(id_o2_yeast) = 0; % Dont allow oxygen secretion

% Use the mattfa solveFBAmodelCplex function to optimize the FBA problem
sol = solveFBAmodelCplex(iMM904);
sol.f % Displat the solution
solutions.yeast_anaerobic = sol.f;


%% Part 3: Phenotypic phase plane

% Ceeate a grid of flux values between 0 and 10
x = 0:1:10; % x for glucose
y = 0:1:10; % y for oxygen
[X,Y] = meshgrid(x,y);

% Allocate an array for the solution
Z = zeros(size(x));

% Loop over the grid
for id_x = 1:length(x)
    for id_y = 1:length(y)
       
       % FIX the glucose and oxygen uptake rates 
       tmodel.lb(id_glc) = -X(id_x,id_y);
       tmodel.ub(id_glc) = -X(id_x,id_y);        
        
       tmodel.lb(id_o2) = -Y(id_x,id_y);
       tmodel.ub(id_o2) = -Y(id_x,id_y);
       
       sol = solveFBAmodelCplex(tmodel);
       
       % Save the solution on a grid
       Z(id_x, id_y) = sol.f;
       
       
    end
end

% Plot the solution
figure_ecoli = surf(X,Y,Z);
xlabel('Glucose uptake [mmol/gDW/hr]')
ylabel('Oxygen uptake [mmol/gDW/hr]')
zlabel('Growth rate [1/hr]')

saveas(figure_ecoli,'E_coli_phase_plane','jpg')


%% Part 4 Phenotypic phase plane with variability

% Create a list of all boundary reactions
% Search for the prefix 'DM'
index_list = [];
for i = 1:length(tmodel.rxns)
    if contains(tmodel.rxns(i), 'DM_') 
        index_list = [index_list , i];
    end 
end

% Ceeate a grid of flux values between 0 and 10 with a step of 5
x = 0:5:10; % x for glucose
y = 0:5:10; % y for oxygen
[X,Y] = meshgrid(x,y);

% Allocate an array for the solution
Z = zeros(size(x));

% Create two 3-dimensional matrices for the variability analysis for all
% the boundary reactions and glucose-oxygen combinations
VA_max = zeros(length(x),length(y),length(index_list));
VA_min = zeros(length(x),length(y),length(index_list));

% Loop over the grid
for id_x = 1:length(x)
    for id_y = 1:length(y)
        
        % Display the position on grid
        id_x
        id_y
        
       % Set the lower bound of biomass to 0 (again see below)
       tmodel.lb(id_biomass) =  0;
       
       % Caculate the optimal growth rate at a single point of 
       % oxygen and glucose uptake 
       tmodel.lb(id_glc) = -X(id_x,id_y);
       tmodel.ub(id_glc) = -X(id_x,id_y);        
        
       tmodel.lb(id_o2) = -Y(id_x,id_y);
       tmodel.ub(id_o2) = -Y(id_x,id_y);
       
       tmodel.c = tmodel.c*0;
       tmodel.c(id_biomass) = 1;
       sol = solveFBAmodelCplex(tmodel);
       
       % store solution in the Z matrix
       Z(id_x, id_y) = sol.f
       
       % Variablity analsysis of the boundary reactions 
       
       % constraint the biomass reaction to the optimal
       tmodel.lb(id_biomass) =  sol.f;
       
       % Loop over the boundary reactions 
       j = 1;
       for sec_id = index_list
           
           % Set the objective to the reaction that we want to perform the
           % variablity analysis 
           tmodel.c = tmodel.c*0;
           tmodel.c(sec_id) = 1;
           
           % Perform maximization of the reaction flux
           sol = solveFBAmodelCplex(tmodel,[],[],[],'max');
           max_value = sol.f;
           % Perform minimization of the reaction flux           
           sol = solveFBAmodelCplex(tmodel,[],[],[],'min');
           min_value = sol.f;
           
           % Store the values of the variability analysis
           VA_max(id_x,id_y,j) = max_value;
           VA_min(id_x,id_y,j) = min_value;
           % > VA_max  = (len_glc, len_o2, len_uptakes (332))
           % > VA_min  = (...) > 0 By products that have to be there 
           %  VA_min = 0  and VA_max > 0 they can be byproducts
           % >>> 3Dim Array 0,1,2 (len_glc, len_o2, len_uptakes (332))
           % FIND unique arrays in Z
           
           j = j + 1;
       end 
       
    end
end

% Based on the previous results we want to identify which metabolites must
% be secreted in each set of oxygen and glucose and which ones can either be
% secreted or not
for id_x = 1:length(x)
    for id_y = 1:length(y)

        % Find the reaction indices for the metabolites that must be secreted
        rxn_idx_forced = find(VA_min(id_x,id_y,:) > 1e-9 );
        % Find the reaction indices for the metabolites that can possibly be secreted
        rxn_idx_can = find((VA_min(id_x,id_y,:) == 0) & (VA_max(id_x,id_y,:) > 0) );        
        % Store the results in cell matrices
        Rxn_forced{id_x,id_y} = tmodel.rxns(rxn_idx_forced);
        Rxn_can{id_x,id_y} = tmodel.rxns(rxn_idx_can);

    end
end

% Save the output (everything that is included in the workspace)
save('./WS_exc1')