%% Example 5: Backtracking of the glycolytic leaks of different species

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

%% Load dynamic model
load('DynModels/Glycolysis_model.mat','DM')

%% Load and prepare GSM Models
%% Yeast %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('GSModels/Yeast/Y7_filled.mat','Y7')
load('DynModels/maps.mat','mapY7')

% disp mapping
mapY7= bktr_mapping(Y7,DM,mapY7);

% set medium conditions (glucose only)
Y7=changeRxnBounds(Y7,'r_1714',-1000,'l');
Y7=changeRxnBounds(Y7,{'r_1549','r_2033','r_1634'},0,'u');


%% Ecoli %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('GSModels/EColi/iJO_filled.mat','iJO1366')
load('DynModels/maps.mat','mapiJO')

% disp mapping
mapiJO= bktr_mapping(iJO1366,DM,mapiJO);

% set medium conditions (glucose only)
iJO1366=changeRxnBounds(iJO1366,'EX_glc(e)',-1000,'l');



%% Mtb   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('GSModels/Mtb/iNJ_filled.mat','iNJ661')
load('DynModels/maps.mat','mapiNJ')

% disp mapping
mapiNJ= bktr_mapping(iNJ661,DM,mapiNJ);

% set medium conditions (glucose only)
iNJ661=changeRxnBounds(iNJ661,{'EX_cit(e)','EX_glu_L(e)','EX_glyc(e)'},-0.001,'l');
iNJ661=changeRxnBounds(iNJ661,'EX_glc(e)',-1000,'l');



%% B.subtilis   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('GSModels/Bsub/iBsu_filled.mat','iBsu1103')
load('DynModels/maps.mat','mapiBsu')

% disp mapping
mapiBsu= bktr_mapping(iBsu1103,DM,mapiBsu);

% set medium conditions (glucose limited)
idxMedium = (findExcRxns(iBsu1103));
iBsu1103.lb(idxMedium)=0;
iBsu1103=changeRxnBounds(iBsu1103,{'EX_cpd00007(e)';'EX_cpd00137(e)';'EX_cpd00001(e)';'EX_cpd10515(e)';'EX_cpd00067(e)';'EX_cpd00058(e)';'EX_cpd00421(e)';'EX_cpd11574(e)';'EX_cpd00531(e)';'EX_cpd11416_c';'EX_cpd10516(e)';'EX_cpd00209(e)';'EX_cpd00971(e)';'EX_cpd00063(e)';'EX_cpd00048(e)';'EX_cpd03424(e)';'EX_cpd00075(e)';'EX_cpd00013(e)';'EX_cpd00244(e)';'EX_cpd00254(e)';'EX_cpd00034(e)';'EX_cpd00149(e)';'EX_cpd00011(e)';'EX_cpd00205(e)';'EX_cpd04097(e)';'EX_cpd00009(e)'},...
    -1000,'l'); %Medium compunds according to Harwood, C. R., and S. M. Cutting. 1990
iBsu1103=changeRxnBounds(iBsu1103,'EX_cpd00027(e)',-1000,'l'); % glucose as sole carbon source

%% Compile all models
models = {Y7 iJO1366  iNJ661  iBsu1103};
map={mapY7 mapiJO  mapiNJ,  mapiBsu};
IdxTargetRxn = [0 8   1025 427];
IdxTargetMet = [ 339   0 0 0 ];
scale = [100 1 50 100 ];
names = {'S.cerevisiae (Y7)','E.coli (iJO1366)','M. tuberculosis (iNJ661)',  'B.Subtilis (iBsu1103)'};


%% Calculate Leaks

% initialize storers
Leaks=zeros(length(mapiJO.DMets), length(models));
Leaks2=Leaks;Leaks_min=Leaks;Leaks_max=Leaks;
AddUptakes=cell(size(models));

for k = 1:length(models)
    
    % Select organism model
    GSM = models{k};
    
    % 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 target
    [GSM, BioRxnIdx]= bktr_Target(GSM,IdxTargetMet(k),IdxTargetRxn(k),scale(k));
    
    % Set up irrev model
    [GSM,irrevGSM,matchRev, map_lump,rev2irrev] = bktr_prep_Lump(GSM,DM,map{k},BioRxnIdx);
    
    % Calculate leaks
    [Leaks(:,k), AddUptakes{k}] = bktr_Run(irrevGSM,GSM,matchRev,map_lump,scale(k));
    
    % Leak FVA (will take some time!) with 5% derivation from the optimum
    % flux dstribution allowed
    [Leaks2(:,k), Leaks_min(:,k), Leaks_max(:,k)]=bktr_LeakFVA(irrevGSM,GSM,matchRev,...
        map_lump,scale(k),rev2irrev,0.95);
    
end

%% Plot results
close all
selected = fliplr([6	7	5	4	8	9	1	2	3	10	11]); % bring in right glycolytic order
bktr_barPlot(Leaks,map_lump.DMets, names)
if any(Leaks2(:)~=0)
    bktr_errorbarPlot(Leaks2(selected,:),Leaks_min(selected,:),Leaks_max(selected,:),map_lump.DMets(selected), names(:))
end





