%% Example 3: Leaks of more example models

%% Initialize

clear all
close all
clc
initCobraToolbox
changeCobraSolver('glpk');
% add all backtracking functions to the searchpath
addpath([cd '/bktr_functions']);

%% Load and prepare GSM and dynamic Models
%% Ecoli Chassagnole %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load GSM model
load('GSModels/EColi/iJO_filled.mat','iJO1366')

% Central carbon metabolism model of chassagnole et al.
DMfile='DynModels/Chassagnole2002/BIOMD0000000051.xml';

DM_SBML = TranslateSBML(DMfile);
DM = convertSBMLToCobra(DM_SBML);
DM = bktr_RefillChebiAndKeggIds(DM,DM_SBML);

% mapping
load('DynModels/Chassagnole2002/map_Chass.mat','mapiJO_Chassangole02')
mapiJO_Chassangole02 = bktr_mapping(iJO1366,DM,mapiJO_Chassangole02);

% set medium conditions
iJO1366=changeRxnBounds(iJO1366,'EX_glc(e)',-1000,'l');  % mmol/l/s

% fix to the original models uptake rates
[~, idx_mup{1}]=ismember({'EX_ac(e)','EX_o2(e)' ,'EX_co2(e)'},iJO1366.rxns);%'EX_glc(e)'
mup{1} = [0  -24.3  24.8 ] /8.7; %mmol/L_ext/h * gDW/L_ext  glucose:-9.2

% Attention: NGAM (ATP costs for maintenence, rxn 'ATPM') are hardcoded as 3.15 mmol/gDW/h in iJO1366


%% Ecoli Singh 2006  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load GSM model
load('GSModels/EColi/iJO_filled.mat','iJO1366')

% TCA cycle on acetate from Singh et al. 2006
DMfile='DynModels/Singh2006/BIOMD0000000221.xml';

DM_SBML = TranslateSBML(DMfile);
DM = convertSBMLToCobra(DM_SBML);
DM = bktr_RefillChebiAndKeggIds(DM,DM_SBML);
DM.metNames= DM.mets;
DM.rxnNames= DM.rxns;

% mapping
load('DynModels/Singh2006/map_Singh.mat','mapiJO_Singh06')
mapiJO_Singh06 = bktr_mapping(iJO1366,DM,mapiJO_Singh06);

% set medium conditions
iJO1366=changeRxnBounds(iJO1366,'EX_glc(e)',-1000,'l');  % mmol/l/s

% fix to the original models uptake rates
[~, idx_mup{2}]=ismember({'EX_glc(e)', 'EX_ac(e)' ,'EX_co2(e)'},iJO1366.rxns);
mup{2} = [0  -7.18  8.76 ] ; %mmol/g/h


%% Yeast Hynne   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% GSM
load('GSModels/Yeast/Y7_filled.mat','Y7');
GSM=Y7;
% DynM
load('DynModels/Hynne_model.mat','DM')
% Map
load('DynModels/maps.mat','mapHynne');
mapHynne = bktr_mapping(Y7,DM,mapHynne);


% set medium conditions
[~, idx_mup{3}]=ismember({'r_1714','r_1631' ,'r_1761', 'r_1808'},GSM.rxns);%glc,aca,eth,glyc
mup{3} = [ -0.82 0.074 0.8 0.073];  %mM/s

% find max growth rate with uptake rates from Hynne et al.
GSM.lb(idx_mup{3})= mup{3};
GSM.ub(idx_mup{3})= mup{3};



%% Yeast van Eunen   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% GSM
load('GSModels/Yeast/Y7_filled.mat','Y7');GSM=Y7;

% Glycolysis from van Eunen et al 2012
DMfile='DynModels/vanEunen2012/MODEL1403250001.xml';

DM_SBML = TranslateSBML(DMfile);
DM = convertSBMLToCobra(DM_SBML);
DM = bktr_RefillChebiAndKeggIds(DM,DM_SBML);
% as no annotation were provided (non curated model in Biomodels), we add the species names from Copasi
DM.metNames= {'Glci';'G6P';'F6P';'F16P';'TRIO';'BPG';'NAD';'NADH';'P3G';'P2G';'PEP';'PYR';'AcAld'};
DM.rxnNames= {'Glt','Hk','Tre1','Tre2','PGI','Pfk','Ald','Gapdh','Gly','Pgk','Gpm','Eno','Pyk','Pdc','Suc','Ace','Adh'};

% mapping
load('DynModels/vanEunen2012/mapY7_vEunen.mat','mapY7_vEunen12')
mapY7_vEunen12 = bktr_mapping(Y7,DM,mapY7_vEunen12);

% set medium conditions
[~, idx_mup{4}]=ismember({'r_1714','r_1634','r_1761','r_1808'},GSM.rxns); %Glucose, Acetate, Ethanol, Glycerol
mup{4} = [ -1.302163 0.1083   2.1743 0.2916] ;  %mM/h


% find max growth rate with uptake rates from van Eunen et al.
GSM.lb(idx_mup{4})= mup{4};
GSM.ub(idx_mup{4})= mup{4};


%% Compile all models
models = {iJO1366 iJO1366 Y7  GSM};
map={mapiJO_Chassangole02 mapiJO_Singh06 mapHynne mapY7_vEunen12};
IdxTargetRxn = [ 8 8 0 1591 ];
IdxTargetMet = [ 0   0 339 0 ];
fix_growth_rate =[0.1 0.11 0  0];%
scale = [1 1 1000  100];
names = {'Chassagnole 2002 (E.coli)','Singh 2006 (E.coli)','Hynne 2001 (Yeast)', 'van Eunen 2012 (Yeast)'};
units={'mmol/(gDW h)', 'mol/(g h)', 'mM/s', 'mM/s'}; %for conversions see Supplementary table 2

%% Calculate Leaks

% initialize storers
Leaks=cell(size(models));
Leaks2=Leaks;Leaks_min=Leaks;Leaks_max=Leaks;
AddUptakes=cell(size(models));

for k = 1:length(models)
    
    % Select organism model
    GSM = models{k};
    
    % Set target
    [GSM, BioRxnIdx]= bktr_Target(GSM,IdxTargetMet(k),IdxTargetRxn(k),scale(k));
    
    if fix_growth_rate(k) ~= 0
        % set growth rate to fixed value (as in the chemostate)
        GSM.lb(BioRxnIdx)=fix_growth_rate(k)/scale(k);
        GSM.ub(BioRxnIdx)=fix_growth_rate(k)/scale(k);
    else
        GSM.lb(BioRxnIdx)=0;
        GSM.ub(BioRxnIdx)=1000;
    end
    
    % measured uptakes
    GSM.lb(idx_mup{k})=mup{k}/scale(k);
    GSM.ub(idx_mup{k})=mup{k}/scale(k);
    
    % adjust reversibility of the fixed fluxes
    GSM.rev([idx_mup{k} BioRxnIdx])=false;
    
    % Display all allowed uptakes
    idxMedium = (findExcRxns(GSM));[a,~]=find(GSM.S(:,idxMedium));
    disp(table(GSM.rxns(idxMedium),printRxnFormula(GSM, GSM.rxns(idxMedium),0,1,1), GSM.lb(idxMedium), GSM.ub(idxMedium),(bktr_nAtoms(GSM,a,'C'))',...
        'VariableNames',{'Rxn','ReactionFormula', 'lb','ub','nC'}))
    
    % Set up irrev model
    [GSM,irrevGSM,matchRev,map_lump,rev2irrev] = bktr_prep_Lump(GSM,DM,map{k},BioRxnIdx);% prepare model
    idup = rev2irrev([idx_mup{k} BioRxnIdx]);
    
    % Calculate leaks
    [Leaks{k}, AddUptakes{k}] = bktr_Run(irrevGSM,GSM,matchRev,map_lump,scale(k));   % Run FBA and Calculate Leaks
    
    subplot(2,2,k), hold on
    bar(1:length(Leaks{k}),Leaks{k})
    set(gca, 'xtick',1:length(Leaks{k}), 'xticklabel', map_lump.DMets,...
        'XTickLabelRotation',45)
    ylabel(['Leak fluxes in ' units{k}])
    title(names(k)), box on, grid on
    
end




