%% 
if strcmp(ModelName, 'yeast8.xml')
    
    GlucoseUptakeReaction = 'r_1714'; % [iaz900-'EX_glc_e_'], [yeast7/'r_1714']
    OxygenUptakeReaction = 'r_1992'; % [iaz900-'EX_o2_e_'], [yeast7/'r_1992']
    ObjectiveRxn = 'r_4041'; % 4041/2111[iaz900-'biomass_wildType'], [yeast7/'r_4041', 'r_2111'] %'BIO0100' %Name of objective reaction in model
    IRON = 'r_1861';
    POTASSIUM = 'r_2020';
    SULPHATE = 'r_2060';
    AMMONIA = 'r_1654';
    PHOSPHATE = 'r_2005';
    SODIUM = 'r_2049';
    WATER = 'r_2100';
    PROTONS = 'r_1832';     

elseif strcmp(ModelName, 'iaz900.xml') == 1
    
    GlucoseUptakeReaction = 'EX_glc_e_'; % [iaz900-'EX_glc_e_'], [yeast7/'r_1714']
    OxygenUptakeReaction = 'EX_o2_e_'; % [iaz900-'EX_o2_e_'], [yeast7/'r_1992']
    ObjectiveRxn = 'biomass_wildType'; % [iaz900-'biomass_wildType'], [yeast7/'r_4041', 'r_2111'] %'BIO0100' %Name of objective reaction in model
    OCozC_name = 'adp[Cytosol]'; % [iaz900-'nad[Cytosol]'], [yeast7/'s_1198[c_03]']
    OCozM_name = 'adp[Mitochondria]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]']
    RCozC_name = 'atp[Cytosol]'; % [iaz900-'nadh[Cytosol]'], [yeast7/'s_1203[c_03]']
    RCozM_name = 'atp[Mitochondria]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    ofset_C = 'pi[Cytosol]';
    ofset_M = 'pi[Mitochondria]';
    IRON = 'EX_fe2_e_';
    POTASSIUM = 'EX_k_e_';
    SULPHATE = 'EX_so4_e_';
    AMMONIA = 'EX_nh4_e_';
    PHOSPHATE = 'EX_pi_e_';
    SODIUM = 'EX_na1_e_';
    WATER = 'EX_h2o_e_';
    PROTONS = 'EX_h_e_';
    
    elseif strcmp(ModelName, 'iJR904.xml') == 1
    
    GlucoseUptakeReaction = 'EX_glc_e_'; % [iaz900-'EX_glc_e_'], [yeast7/'r_1714']
    OxygenUptakeReaction = 'EX_o2_e_'; % [iaz900-'EX_o2_e_'], [yeast7/'r_1992']
    ObjectiveRxn = 'BiomassEcoli'; % [iaz900-'biomass_wildType'], [yeast7/'r_4041', 'r_2111'] %'BIO0100' %Name of objective reaction in model
    OCozC_name = 'nad[Cytosol]'; % [iaz900-'nad[Cytosol]'], [yeast7/'s_1198[c_03]']
    %OCozM_name = 'nad[Mitochondria]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]']
    RCozC_name = 'nadh[Cytosol]'; % [iaz900-'nadh[Cytosol]'], [yeast7/'s_1203[c_03]']
    %RCozM_name = 'nadh[Mitochondria]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    ofset_C = '';
    ofset_M = '';
    IRON = 'EX_fe2_e_';
    POTASSIUM = 'EX_k_e_';
    SULPHATE = 'EX_so4_e_';
    AMMONIA = 'EX_nh4_e_';
    PHOSPHATE = 'EX_pi_e_';
    SODIUM = 'EX_na1_e_';
    WATER = 'EX_h2o_e_';
    PROTONS = 'EX_h_e_';
    
elseif strcmp(ModelName, 'iCEL1273.xml') == 1
    
    GlucoseUptakeReaction = 'EX00031'; % [iaz900-'EX_glc_e_'], [yeast7/'r_1714']
    OxygenUptakeReaction = 'EX00007'; % [iaz900-'EX_o2_e_'], [yeast7/'r_1992']
    ObjectiveRxn = 'BIO0100'; % [iaz900-'biomass_wildType'], [yeast7/'r_4041', 'r_2111'] %'BIO0100' %Name of objective reaction in model
    OCozC_name = 'C00003[Cytosol]'; % [iaz900-'nad[Cytosol]'], [yeast7/'s_1198[c_03]']
    OCozM_name = 'M00003[Mitochondria]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]']
    RCozC_name = 'C00004[Cytosol]'; % [iaz900-'nadh[Cytosol]'], [yeast7/'s_1203[c_03]']
    RCozM_name = 'M00004[Mitochondria]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    ofset_C = '';
    ofset_M = '';
    PHOSPHATE = 'EX00009';
    WATER = 'EX00001';
    PROTONS = 'EX00080';
    BacterialUptakeReaction = 'EXC0050';
    
elseif strcmp(ModelName, 'iND750.xml') == 1
    
    GlucoseUptakeReaction = 'EX_glc_e_'; % [iaz900-'EX_glc_e_'], [yeast7/'r_1714']
    OxygenUptakeReaction = 'EX_o2_e_'; % [iaz900-'EX_o2_e_'], [yeast7/'r_1992']
    ObjectiveRxn = 'biomass_SC4_bal'; % [iaz900-'biomass_wildType'], [yeast7/'r_4041', 'r_2111'] %'BIO0100' %Name of objective reaction in model
    OCozC_name = 'nad[Cytosol]'; % [iaz900-'nad[Cytosol]'], [yeast7/'s_1198[c_03]']
    OCozM_name = 'nad[Mitochondria]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]']
    RCozC_name = 'nadh[Cytosol]'; % [iaz900-'nadh[Cytosol]'], [yeast7/'s_1203[c_03]']
    RCozM_name = 'nadh[Mitochondria]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    ofset_C = '';
    ofset_M = '';
    IRON = 'EX_fe2_e_';
    POTASSIUM = 'EX_k_e_';
    SULPHATE = 'EX_so4_e_';
    AMMONIA = 'EX_nh4_e_';
    PHOSPHATE = 'EX_pi_e_';
    SODIUM = 'EX_na1_e_';
    WATER = 'EX_h2o_e_';
    PROTONS = 'EX_h_e_';
elseif strcmp(Organism, 'H_sapiens') == 1    
    GlucoseUptakeReaction = 'EX_glc(e)'; % 'EX_glc(e)' glucose, 'EX_sucr(e)' sucrose
    OxygenUptakeReaction = 'EX_o2(e)'; % [iaz900-'EX_o2_e_'], [yeast7/'r_1992']
    ObjectiveRxn = 'biomass_reaction'; %Name of objective reaction in model
    %  'DM_atp_c_' is atp maintenance; biomass_reaction is biomass  
    OCozC_name = 'nad[c]'; % [iaz900-'nad[Cytosol]'], [yeast7/'s_1198[c_03]']
    RCozC_name = 'nadh[c]'; % [iaz900-'nadh[Cytosol]'], [yeast7/'s_1203[c_03]']
    OCozM_name = 'nad[m]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]'
    RCozM_name = 'nadh[m]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    OCozN_name = 'nad[n]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]'
    RCozN_name = 'nadh[n]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    OCozE_name = 'nad[e]'; % [iaz900-'nad[Mitochondria]'], [yeast7/'s_1200[c_11]'
    RCozE_name = 'nadh[e]'; % [iaz900-'nadh[Mitochondria]'], [yeast7/'s_1205[c_11]']
    ofset_C = '';
    ofset_M = '';
    IRON = 'EX_fe3(e)';
    SELENIUM = 'EX_sel(e)';
    IODINE = 'EX_i(e)';
%     COPPER = ;
%     MAGNESIUM = ;
%     MANGANESE = ;
    ZINC = 'EX_HC02172(e)';
    POTASSIUM = 'EX_k(e)';
    SODIUM = 'EX_na1(e)';
    PHOSPHATE = 'EX_pi(e)';
    CALCIUM = 'EX_ca2(e)';
    CHLORIDE = 'EX_cl(e)';
%    'EX_fe2(u)'     
    SULPHITE = 'EX_so3(e)';
    NITRITE = 'EX_no2(e)';
    AMMONIA = 'EX_nh4(e)';
    WATER = 'EX_h2o(e)';
    PROTONS = 'EX_h(e)';
    HYDROXIDE = 'EX_oh1(e)';
    BICARBONATE = 'EX_hco3(e)';
    CO2 = 'EX_co2(e)';
    EA1_histidine = 'EX_his_L(e)';
    EA2_isoleucine = 'EX_ile_L(e)';
    EA3_leucine = 'EX_leu_L(e)';
    EA4_lysine = 'EX_lys_L(e)';
    EA5_methionine = 'EX_met_L(e)';
    EA6_phenylalanine = 'EX_phe_L(e)';
    EA7_threonine = 'EX_thr_L(e)';
    EA8_tryptophan = 'EX_trp_L(e)';
    EA9_valine = 'EX_val_L(e)';
    VA_retinol = 'EX_11_cis_retfa(e)'; %'EX_retinol(e)'
    VB1_thyamine = 'EX_thm(e)';
    VB2_riboflavin = 'EX_ribflv(e)';
    VB3_niacins = 'DM_ncam';
    VB3_niacins_NAM = 'EX_ncam(e)';
    VB3_niacins_NA = 'EX_nac(e)';
    VB3_niacins_NMNA = 'EX_1mncam(e)';
%     VB3_niacins_NMN = 'DM_ncam';
%     VB3_niacins_NR = 'DM_ncam';
    VB5_panthotenate = 'DM_pnto_R';
    VB6_pyridoxine = 'EX_avite1(e)';
    VB7_biotin = 'DM_btn';
    VB8_myoinositol = 'EX_inost(e)';
    VB8_inosine = 'EX_ins(e)';
    VB9_folate = 'DM_fol';
    VB12_cobalamine = 'EX_adpcbl(e)';
    VC_ascorbate = 'EX_ascb_L(e)';
    proVD_colecalciferol = 'EX_vitd2(e)';
    proVD_D3 = 'EX_vitd3(e)';
    VE_tocopherol = 'EX_avite1(e)';
    VK_naphtoquinones_phyQ = 'EX_phyQ(e)';
    
else
end

%% 

model = readCbModel();

    Exs_loc = findExcRxns(model);
    Exs = model.rxns(Exs_loc);
    model = changeRxnBounds(model, Exs, 0, 'l');
    if strcmp(Organism, 'C_elegans') == 1
            model = changeRxnBounds(model, WATER, -1000, 'l');
            model = changeRxnBounds(model, PROTONS, -1000, 'l');
            model = changeRxnBounds(model, PHOSPHATE, -1000, 'l');
            model = changeRxnBounds(model, BacterialUptakeReaction, -0.32, 'l');
            
        elseif strcmp(Organism, 'H_sapiens') == 1
            model = changeRxnBounds(model, IRON, -1000, 'l');
            model = changeRxnBounds(model, SELENIUM, -1000, 'l');
            model = changeRxnBounds(model, IODINE, -1000, 'l');
           
            model = changeRxnBounds(model, ZINC, -1000, 'l');

            model = changeRxnBounds(model, POTASSIUM, -1000, 'l');
            model = changeRxnBounds(model, SODIUM, -1000, 'l');
            model = changeRxnBounds(model, PHOSPHATE, -1000, 'l');
            model = changeRxnBounds(model, CALCIUM, -1000, 'l');
            model = changeRxnBounds(model, CHLORIDE, -1000, 'l'); 
            %model = changeRxnBounds(model, SULPHITE, -1000, 'l');
            %model = changeRxnBounds(model, NITRITE, -1000, 'l');
            %model = changeRxnBounds(model, AMMONIA, -1000, 'l');
            model = changeRxnBounds(model, WATER, -1000, 'l');
            model = changeRxnBounds(model, PROTONS, -1000, 'l');
            model = changeRxnBounds(model, HYDROXIDE, -1000, 'l');
            model = changeRxnBounds(model, BICARBONATE, -1000, 'l');
            model = changeRxnBounds(model, CO2, -1000, 'l');
            model = changeRxnBounds(model, EA1_histidine, -1000, 'l');
            model = changeRxnBounds(model, EA2_isoleucine, -1000, 'l');
            model = changeRxnBounds(model, EA3_leucine, -1000, 'l');
            model = changeRxnBounds(model, EA4_lysine, -1000, 'l');
            model = changeRxnBounds(model, EA5_methionine, -1000, 'l');
            model = changeRxnBounds(model, EA6_phenylalanine, -1000, 'l');
            model = changeRxnBounds(model, EA7_threonine, -1000, 'l');
            model = changeRxnBounds(model, EA8_tryptophan, -1000, 'l');
            model = changeRxnBounds(model, EA9_valine, -1000, 'l');
            model = changeRxnBounds(model, VA_retinol, -1000, 'l');
            model = changeRxnBounds(model, VB1_thyamine, -1000, 'l');
            model = changeRxnBounds(model, VB2_riboflavin, -1000, 'l');
            model = changeRxnBounds(model, VB3_niacins, -1000, 'l');
            model = changeRxnBounds(model, VB3_niacins_NAM, 0, 'l');
            model = changeRxnBounds(model, VB3_niacins_NA, 0, 'l');
            model = changeRxnBounds(model, VB3_niacins_NMNA, 0, 'l');
            model = changeRxnBounds(model, VB5_panthotenate, -1000, 'l');
            model = changeRxnBounds(model, VB6_pyridoxine, -1000, 'l');
            model = changeRxnBounds(model, VB7_biotin, -1000, 'l');
            model = changeRxnBounds(model, VB8_myoinositol, -1000, 'l');
            model = changeRxnBounds(model, VB8_inosine, -1000, 'l');
            model = changeRxnBounds(model, VB9_folate, -1000, 'l');
            model = changeRxnBounds(model, VB12_cobalamine, -1000, 'l');
            model = changeRxnBounds(model, VC_ascorbate, -1000, 'l');
            model = changeRxnBounds(model, proVD_colecalciferol, -1000, 'l');
            model = changeRxnBounds(model, proVD_D3, -1000, 'l');
            model = changeRxnBounds(model, VE_tocopherol, -1000, 'l');
            model = changeRxnBounds(model, VK_naphtoquinones_phyQ, -1000, 'l');

    else
        model = changeRxnBounds(model, IRON, -1000, 'l');
        model = changeRxnBounds(model, POTASSIUM, -1000, 'l');
        model = changeRxnBounds(model, SODIUM, -1000, 'l');
        model = changeRxnBounds(model, WATER, -1000, 'l');
        model = changeRxnBounds(model, PROTONS, -1000, 'l');
        model = changeRxnBounds(model, SULPHATE, -1000, 'l');
        model = changeRxnBounds(model, AMMONIA, -1000, 'l');
        model = changeRxnBounds(model, PHOSPHATE, -1000, 'l');
    end


model = buildRxnGeneMat(model);
model = changeObjective(model, ObjectiveRxn);

fba = optimizeCbModel(model, 'max', 'one')
if strcmp(Medium, 'Unrestrained') == 1      
    model = changeRxnBounds(model, Exs, -1000, 'l');
elseif strcmp(ModelName, 'iaz900.xml') == 1 && strcmp(Medium, 'Rich') == 1

 ObjectiveRxn_In = find(strcmp(model.rxns, ObjectiveRxn) == 1);
 Cytosol = '[Cytosol]';
 Extra = '[Extra_organism]';
 Exs_In = 486:652;
 Blocks = findMetsFromRxns(model, ObjectiveRxn_In);
    for el = 1: length(Blocks)
            Blocks(el) = erase(Blocks(el), Cytosol);
            Blocks(el) = strcat(Blocks(el), Extra);
    end

   Sugars{1, 1} = 'glc[Extra_organism]';
   Sugars{2, 1} = 'g6p[Extra_organism]';
        Sugars{3, 1} = 'man6p[Extra_organism]';
        Sugars{4, 1} = 'mnl[Extra_organism]';
        Sugars{5, 1} = 'lcts[Extra_organism]';
        Sugars{6, 1} = 'sucr[Extra_organism]';
        Sugars{7, 1} = 'malt[Extra_organism]';
        Sugars{8, 1} = 'fru[Extra_organism]';
        Sugars{9, 1} = 'sbt_D[Extra_organism]';
        Sugars{10, 1} = 'man[Extra_organism]';
        Sugars{11, 1} = 'glcn[Extra_organism]';
        Sugars{12, 1} = 'arab_L[Extra_organism]';
        Sugars{13, 1} = 'rib_D[Extra_organism]';
        Sugars{14, 1} = 'xyl_D[Extra_organism]';
        Sugars{15, 1} = 'gal[Extra_organism]';
        
        Basics{1, 1} = 'pi[Extra_organism]';
        Basics{2, 1} = 'nh4[Extra_organism]';
        Basics{3, 1} = 'na1[Extra_organism]';
        Basics{4, 1} = 'k[Extra_organism]';
        Basics{5, 1} = 'h[Extra_organism]';
        Basics{6, 1} = 'h2o[Extra_organism]';
        Basics{7, 1} = 'fe2[Extra_organism]';
        Basics{8, 1} = 'o2[Extra_organism]';
        Basics{9, 1} = 'co2[Extra_organism]';

   Exs_nongrIn = zeros(size(Exs));
    for el = 1 : length(Exs)
        Nutrients = findMetsFromRxns(model, Exs_In(el));
        cont = zeros(1, length(Nutrients));
        for ol = 1 : length(Nutrients)
            if any(strcmp(Blocks, Nutrients(ol))==1) || any(strcmp(Sugars, Nutrients(ol))==1) || any(strcmp(Basics, Nutrients(ol))==1)
                cont(ol) = 1;
            end
        end
        if sum(cont) == 0
            Exs_nongrIn(el) = Exs_In(el);
        end
    end
    
    Exs_nongrIn = Exs_nongrIn(Exs_nongrIn ~= 0);
    Exs_nongr = model.rxns(Exs_nongrIn);
    
    model = changeRxnBounds(model, Exs_nongr, -1000, 'l');
end
       

if strcmp(Optimization, 'FBA') == 1
    model = changeRxnBounds(model, GlucoseUptakeReaction, GlucoseUptake, 'l');
    model = changeRxnBounds(model, OxygenUptakeReaction, OxygenUptake, 'l');
end

genome = model.genes;
Ng = length(genome);
%%
if strcmp(Optimization, 'CAFBA') == 1
    if strcmp(ModelName, 'iaz900.xml') == 1
        % Remove default bounds for glucose and oxygen
        ex_glc_r=find(strcmp(model.rxns,'EX_glc_e_'));
        ex_o2_r=find(strcmp(model.rxns,'EX_o2_e_'));
        model.lb(ex_glc_r)=GlucoseUptake; model.ub(ex_glc_r)=1000;
        model.lb(ex_o2_r)=OxygenUptake; model.ub(ex_o2_r)=1000;
        % Knocking out the GLCDe flux
        glcd_r=find(strcmp(model.rxns,'GLCDe'));
        model.lb(glcd_r)=0;
        model.ub(glcd_r)=0;
        model = addProteinGroupsToModel_iaz900(model);
%         
%         % We now set the offsets for the different groups.
%         model.protGroup(1).phi0 = 0;
%         model.protGroup(2).phi0 = 0;
%         model.protGroup(3).phi0 = 0.08;
%         model.protGroup(4).phi0 = 0.5471;
%         
%         % Now we have to set the weights for the different groups.
%         % (note that the function creates the field model.weights if needed)
%         model=setWeights(model,1,0.015);
%         model=setWeights(model,2,0.00083);
%         model=setWeights(model,3,0.35);
%         model=setWeights(model,4,0);

        % Clear useless variables
        clear('ex_glc_r','ex_o2_r');
       
    elseif strcmp(ModelName, 'yeast7.xml') == 1
                % Remove default bounds for glucose and oxygen
                ex_glc_r=find(strcmp(model.rxns,'r_1714'));
                ex_o2_r=find(strcmp(model.rxns,'r_1992'));
                model.lb(ex_glc_r)=GlucoseUptake; model.ub(ex_glc_r)=1000;
                model.lb(ex_o2_r)=OxygenUptake; model.ub(ex_o2_r)=1000;
                % Knocking out the GLCDe flux
                glcd_r=find(strcmp(model.rxns,'GLCDe'));
                model.lb(glcd_r)=0;
                model.ub(glcd_r)=0;
                model = addProteinGroupsToModel_yeast7(model);
                %
                %         % We now set the offsets for the different groups.
                %         model.protGroup(1).phi0 = 0;
                %         model.protGroup(2).phi0 = 0;
                %         model.protGroup(3).phi0 = 0.08;
                %         model.protGroup(4).phi0 = 0.5471;
                %
                %         % Now we have to set the weights for the different groups.
                %         % (note that the function creates the field model.weights if needed)
                %         model=setWeights(model,1,0.015);
                %         model=setWeights(model,2,0.00083);
                %         model=setWeights(model,3,0.35);
                %         model=setWeights(model,4,0);
                
                % Clear useless variables
                clear('ex_glc_r','ex_o2_r');
    end
    model = setYeastCtAlocationParameters(model, GlucoseUptake, OxygenUptake);
    
    clear('ex_glc_r','ex_o2_r');
end
%%
Name_Imbalance1 = strcat('Imbalance_', PMETag, '_', ArrangementName_1);
if strcmp(MitFlag, 'ON') == 1
Name_Imbalance2 = strcat('Imbalance_', PMETag, '_', ArrangementName_2);
end

model_redox = addReaction(model, Name_Imbalance1, {RCozC_name, OCozC_name, ofset_C}, [-1 1 1]);

if strcmp(MitFlag, 'ON') == 1
model_redox = addReaction(model_redox, Name_Imbalance2, {RCozM_name, OCozM_name, ofset_M}, [-1 1 1]);
end

clc