%% Example 6: Modification of model topologies


%% Initialize
clear all, close all, clc
initCobraToolbox
changeCobraSolver('glpk');

% add all backtracking functions to the searchpath
addpath([cd '/bktr_functions']);


%% Load GSM and initial DDM and initialize storers for the modified models
DMall= cell(1,5);
map_all = cell(1,5);


load('DynModels/Glycolysis_model.mat','DM')
load('GSModels/Yeast/Y7_filled.mat','Y7')


%% Compile model variants
%% Dynamic model 1: original dummy glycolysis model
load('DynModels/maps.mat','mapY7')
map= bktr_mapping(Y7,DM,mapY7);

DMall{1} = DM;
map_all{1} = map;


%% Dynamic model 2&3: Splitting in upper and lower glycolysis

idx_rxn{2} =[3 4 6 1 2 11]; idx_rxn{3}=[7 8 9 10];
idx_met{2} =[4 5 6 7 8 9 12:13]; idx_met{3} =[1 2 3 10 11 12:13];

for k=2:3
    %initialize
    DMall{k} = DM;
    map_all{k} = mapY7;
        
    idx_rxns=idx_rxn{k};
    idx_mets=idx_met{k};
    
    % fill map
    map_all{k}.DRxns = map.DRxns(idx_rxns);
    map_all{k}.DMets = map.DMets(idx_mets);
    map_all{k}.GSMRxns=map.GSMRxns(idx_rxns);
    map_all{k}.GSMMets=map.GSMMets(idx_mets);

    % fill model
    DMall{k}.mets = DM.mets(idx_mets);
    DMall{k}.rxns = DM.rxns(idx_rxns);
    DMall{k}.metNames = map.DMets(idx_mets);
    DMall{k}.rxnNames = map.DRxns(idx_rxns);
    DMall{k}.S=DM.S(idx_mets,idx_rxns);
    DMall{k}.rev = DM.rev(idx_rxns);
    DMall{k}.metKEGGID= DM.metKEGGID(idx_mets);
    DMall{k}.metChEBIID=cell(size(DM.mets(idx_mets)));
    DMall{k}.rxnKEGGID= DM.rxnKEGGID(idx_rxns);
    DMall{k}.metFormula=DM.metFormula(idx_mets);

    % show mapping
    map_all{k}=bktr_mapping(Y7,DMall{k},map_all{k});
     
end

%% Dynamic model 4: Simplified by lumping reactions
k=4;

% fill map
map_all{k}.DRxns = {'upper glycolysis';'fructose-bisphosphate aldolase';'triose-phosphate isomerase';'lower glycolysis start';'lower glycolysis pyr production'};
map_all{k}.DMets = {'D-glucose';'D-fructose 1,6-bisphosphate';'glyceraldehyde 3-phosphate';'dihydroxyacetone phosphate';'phospoglycerate  pool';'pyruvate';'ADP';'ATP';'NAD';'NADH'};
YRxns={{'r_0534','r_0467','r_0886','r_0449'};'r_0450';'r_1054';{'r_0486','r_0892'};{'r_0366','r_0962'}};
map_all{k}.GSMRxns=cell(size(YRxns));
for i=1:length(YRxns)
    map_all{k}.GSMRxns{i}=findRxnIDs(Y7,YRxns{i});
end
YMets={'s_0563';'s_0555';'s_0764';'s_0629';{'s_0188','s_0260'};'s_1399';'s_0394';'s_0434';'s_1198';'s_1203'};
map_all{k}.GSMMets=cell(size(YRxns));
for i=1:length(YMets)
    map_all{k}.GSMMets{i}=findMetIDs(Y7,YMets{i});
end

% fill model
DMall{k}.mets = {'glc-D[c]';'f6p[c]';'g3p[c]';'dhap[c]';'xpg[c]';'pyr[c]';'adp[c]';'atp[c]';'nad[c]';'nadh[c]'};
DMall{k}.rxns = {'LR1';'LR2';'TPI';'LR3';'LR4'};
DMall{k}.metNames = map_all{k}.DMets;
DMall{k}.rxnNames = map_all{k}.DRxns;
DMall{k}.S=sparse([-1,0,0,0,0;1,-1,0,0,0;0,1,1,-1,0;0,1,-1,0,0;0,0,0,1,-1;0,0,0,0,1;1,1,0,-1,-1;-1,-1,0,1,1;0,0,0,-1,0;0,0,0,1,0]);
DMall{k}.rev = [1;1;1;1;1];
DMall{k}.metKEGGID= {'C00031';'C00085';'C00661';'C00111';{'C00631','C00197'};'C00022';'C00008';'C00002';'C00003';'C00004'};
DMall{k}.metChEBIID=cell(size(DMall{k}.mets));
DMall{k}.rxnKEGGID= {{'R00332','R02003','R01377'};'R00519';'R02016';{'R00115','R00756'};{'R01009',[]}};
DMall{k}.metFormula={'C6H12O6';'C6H13O9P';'C3H7O6P';'C3H7O6P';{'C3H7O7P','C3H7O7P'};'C3H4O3';'C10H15N5O10P2';'C10H16N5O13P3';'C21H28N7O14P2';'C21H29N7O14P2'};

% show mapping
map_all{k}=bktr_mapping(Y7,DMall{k},map_all{k});


%% Dynamic model 5: Extended by pentose phosphate pathway
k=5;

% additional reactions and metabolites
add_mets={'s_0340','s_0335','s_0551','s_0577','s_1408','s_0581','s_1427','s_1207','s_1212'};
add_rxns={'r_0091','r_0466','r_0889','r_0982','r_0984','r_1048','r_1049','r_1050'};
idx_add_mets = findMetIDs(Y7,add_mets);
idx_add_rxns = findRxnIDs(Y7,add_rxns);

% fill map
map_all{k}.GSMRxns=[map.GSMRxns num2cell(idx_add_rxns)];
map_all{k}.GSMMets=[map.GSMMets num2cell(idx_add_mets)];
map_all{k}.DRxns =[ map.DRxns ; Y7.rxnNames(idx_add_rxns)];
map_all{k}.DMets =[ map.DMets ; Y7.metNames(idx_add_mets)];

% fill model
DMall{k}.mets = [DM.mets  Y7.mets(idx_add_mets)' ];
DMall{k}.rxns = [DM.rxns  Y7.rxns(idx_add_rxns)' ];
DMall{k}.metNames = Y7.metNames([map.GSMMets{:} idx_add_mets]);
splitnames=cellfun(@(x) strsplit(x, ' ['), DMall{k}.metNames, 'uniformOutput', false);
splitnames=[splitnames{:}];
DMall{k}.metNames=splitnames(1:2:end)';
DMall{k}.rxnNames = Y7.rxnNames([map.GSMRxns{:} idx_add_rxns]);
DMall{k}.S=Y7.S([map_all{k}.GSMMets{:}],[map_all{k}.GSMRxns{:}]);
DMall{k}.rev = [DM.rev ; Y7.rev(idx_add_rxns)];
DMall{k}.metKEGGID= [DM.metKEGGID ;Y7.metKEGGID(idx_add_mets)];
DMall{k}.metChEBIID=cell(size(DMall{k}.mets'));
DMall{k}.rxnKEGGID= [DM.rxnKEGGID ; cell(size(idx_add_rxns))'];
DMall{k}.metFormula=[DM.metFormula ; Y7.metFormula(idx_add_mets)];

% show mapping
map_all{k}=bktr_mapping(Y7,DMall{k},map_all{k});


%% set medium conditions for the GSM (glucose only)
Y7=changeRxnBounds(Y7,'r_1714',-1000,'l');

%% Define Targtes
IdxTargetRxn = 0;
IdxTargetMet = 339 ;
scale = 1000; 


%% Run Backtracking with different dynamic models
for k=1:length(DMall)
    
    % Select dynamic model
    map = map_all{k};
    DM = DMall{k};
    
    % set target
    [GSM, BioRxnIdx]= bktr_Target(Y7,IdxTargetMet,IdxTargetRxn,scale);
    
    % Prepare and run Backtracking
    [GSM,irrevGSM,matchRev, map_lump,rev2irrev] = bktr_prep_Lump(GSM,DM,map,BioRxnIdx);% prepare model
    [Leaks{k}, AddUptakes{k},waste{k}] = bktr_Run(irrevGSM,GSM,matchRev,map_lump,scale);   % Run FBA and Calculate Leaks
    
    
end


%% Display results
CollectLeaks=NaN(length(DMall{end}.metNames),length(DMall));
for k= 1:length(DMall)
    [a,b]=ismember(DMall{k}.metNames,DMall{end}.metNames); 
    CollectLeaks(b(a),k)=Leaks{k}(a);
end
results = table(DMall{end}.metNames, CollectLeaks,...
    'VariableNames',{'Metabolites','Leaks'});



