setColorToMapRxns

 USAGE:
 [modifiedMap, xmlMap, model2] =
   setColorToMapRxns (mapFileName, model, v1, rxnsFluxTask, rxnLineWidth, rxnLineColor, molFillColor)

 INPUTS:
 mapFileName           relative path to XML map (e.g. 'yeastMap1.0.xml')
 model                 model structure
 v1                    vector of fluxes
 rxnsFluxTask          1 = Flux
                       2 = Flux essentiality; map reactions causing any
                           effect on growth (i.e. obj function)
                       3 = Flux essentiality; map reactions causing no
                           growth(i.e. obj function) when deleted, i.e.
                           essential reactions
                       4 = Flux essentiality; map reactions with not any
                           effect on growth (i.e. obj function) when
                           deleted
                       5 = Omics
 ADDITIONAL INPUTS:
 rxnLineWidth          reaction line width; default = 5
 rxnLineColor          reaction line color; default = Pathways; write
                       'Subsystems' to color five subsystems only (central,
                       AA, lipids, cofactors, nuclotides and energy).
                       'Subsystems' option is the deault option at task 5
 molFillColor          molecule fill color; default = true (i.e. white)
                       molecule fill color; true (i.e. like CellDesigner default)
 OUTPUTS:
 modifiedMap           map with MATLAB structure, and modified reactions' line color and width
 xmlMap                xml structure of the original map obtained from mapFileName
 model2                model with modified subsystems

0001 function [modifiedMap, xmlMap, model2] = setColorToMapRxns (mapFileName, model, v1, rxnsFluxTask, rxnLineWidth, rxnLineColor, molFillColor)
0002 % USAGE:
0003 % [modifiedMap, xmlMap, model2] =
0004 %   setColorToMapRxns (mapFileName, model, v1, rxnsFluxTask, rxnLineWidth, rxnLineColor, molFillColor)
0005 %
0006 % INPUTS:
0007 % mapFileName           relative path to XML map (e.g. 'yeastMap1.0.xml')
0008 % model                 model structure
0009 % v1                    vector of fluxes
0010 % rxnsFluxTask          1 = Flux
0011 %                       2 = Flux essentiality; map reactions causing any
0012 %                           effect on growth (i.e. obj function)
0013 %                       3 = Flux essentiality; map reactions causing no
0014 %                           growth(i.e. obj function) when deleted, i.e.
0015 %                           essential reactions
0016 %                       4 = Flux essentiality; map reactions with not any
0017 %                           effect on growth (i.e. obj function) when
0018 %                           deleted
0019 %                       5 = Omics
0020 % ADDITIONAL INPUTS:
0021 % rxnLineWidth          reaction line width; default = 5
0022 % rxnLineColor          reaction line color; default = Pathways; write
0023 %                       'Subsystems' to color five subsystems only (central,
0024 %                       AA, lipids, cofactors, nuclotides and energy).
0025 %                       'Subsystems' option is the deault option at task 5
0026 % molFillColor          molecule fill color; default = true (i.e. white)
0027 %                       molecule fill color; true (i.e. like CellDesigner default)
0028 % OUTPUTS:
0029 % modifiedMap           map with MATLAB structure, and modified reactions' line color and width
0030 % xmlMap                xml structure of the original map obtained from mapFileName
0031 % model2                model with modified subsystems
0032 
0033 if ~(exist('transformXML2Map')==2)
0034     error('COBRA Toolbox is required, go to https://opencobra.github.io/cobratoolbox/')
0035 end
0036 
0037 if nargin < 5 || rxnLineWidth == 0
0038     rxnLineWidth = 0;
0039     rxnLineColor = '';
0040     molFillColor = true;
0041 elseif nargin < 6
0042     rxnLineColor = '';
0043     molFillColor = true;
0044 elseif nargin < 7
0045     molFillColor = true;
0046 end
0047 if rxnsFluxTask == 5
0048     rxnLineColor = 'Subsystems';
0049 end
0050 
0051 [xmlMap, mapMap] = transformXML2Map(mapFileName);
0052 if molFillColor == true
0053     modifiedMap = unifyMetabolicMapCD(mapMap); % molColor whrite and rxnColor = ligthgray
0054 elseif molFillColor == false
0055     modifiedMap = mapMap; % like CellDesigner default
0056 end
0057 
0058 %%%%%%%% Gettin reactions line width with normalized fluxes %%%%%%%%%%%%%
0059 model2 = model; % to not overwrite model
0060 flxs(:, 1) = v1; % to not overwrite v1
0061 absFlxs = abs(flxs);
0062 rxnLineWidthInMap = absFlxs / max(absFlxs); %normalize fluxes line width
0063 if rxnsFluxTask == 1 && rxnLineWidth == 0
0064     rxnLineWidthInMap(rxnLineWidthInMap > 0.1 & rxnLineWidthInMap <=  1) = 50;
0065     rxnLineWidthInMap(rxnLineWidthInMap > 0.01 & rxnLineWidthInMap <= 0.1) = 25;
0066     rxnLineWidthInMap(rxnLineWidthInMap > 0.001 & rxnLineWidthInMap <= 0.01) = 16;
0067     rxnLineWidthInMap(rxnLineWidthInMap > 0.0001 & rxnLineWidthInMap <= 0.001) = 10;
0068     rxnLineWidthInMap(rxnLineWidthInMap > 0.00001 & rxnLineWidthInMap <= 0.0001) = 6;
0069 elseif rxnsFluxTask == 2 && rxnLineWidth == 0
0070     rxnLineWidthInMap(rxnLineWidthInMap > 0.90 & rxnLineWidthInMap <=  0.95) =  6;
0071     rxnLineWidthInMap(rxnLineWidthInMap > 0.80 & rxnLineWidthInMap <= 0.90) =  8;
0072     rxnLineWidthInMap(rxnLineWidthInMap > 0.50 & rxnLineWidthInMap <= 0.80) = 11;
0073     rxnLineWidthInMap(rxnLineWidthInMap > 0.25 & rxnLineWidthInMap <= 0.50) = 15;
0074     rxnLineWidthInMap(rxnLineWidthInMap >= 0.00 & rxnLineWidthInMap <= 0.25) = 20;
0075 elseif rxnsFluxTask == 1 && rxnLineWidth ~= 0.0
0076     rxnLineWidthInMap(rxnLineWidthInMap == 1.0) = rxnLineWidth;
0077 elseif rxnsFluxTask == 2 && rxnLineWidth ~= 0.0
0078     rxnLineWidthInMap(rxnLineWidthInMap < 0.95) = rxnLineWidth;
0079 elseif rxnsFluxTask == 3 && rxnLineWidth == 0.0
0080     rxnLineWidthInMap(rxnLineWidthInMap == 0.0) = 5;
0081 elseif rxnsFluxTask == 3 && rxnLineWidth ~= 0.0
0082     rxnLineWidthInMap(rxnLineWidthInMap == 0.0) = rxnLineWidth;
0083 elseif rxnsFluxTask == 4 && rxnLineWidth == 0.0
0084     rxnLineWidthInMap(rxnLineWidthInMap > 0.95 & rxnLineWidthInMap <= 1.0) = 5.0;
0085 elseif rxnsFluxTask == 4 && rxnLineWidth ~= 0.0
0086     rxnLineWidthInMap(rxnLineWidthInMap > 0.95 & rxnLineWidthInMap <= 1.0) = rxnLineWidth;
0087 elseif rxnsFluxTask == 5 && rxnLineWidth == 0.0
0088     cutValue = 2.0; % transcriptomics minimum cut value of log change
0089     rxnLineWidthInMap = absFlxs;
0090     rxnLineWidthInMap(rxnLineWidthInMap > 5.0) = 30;
0091     rxnLineWidthInMap(rxnLineWidthInMap > 4.0 & rxnLineWidthInMap <= 5.0) = 25;
0092     rxnLineWidthInMap(rxnLineWidthInMap > 3.0 & rxnLineWidthInMap <= 4.0) = 20;
0093     rxnLineWidthInMap(rxnLineWidthInMap > 2.0 & rxnLineWidthInMap <= 3.0) = 15;
0094     rxnLineWidthInMap(rxnLineWidthInMap > 1.0 & rxnLineWidthInMap <= 2.0) = 10;
0095     rxnLineWidthInMap(rxnLineWidthInMap >= 0.5 & rxnLineWidthInMap <= 1.0) =  5;
0096 end
0097 for i = 1:length(modifiedMap.rxnName)
0098     a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0099     if isempty(a)
0100         modifiedMap.rxnWidth{i} = 1;
0101     else
0102         modifiedMap.rxnWidth{i} = rxnLineWidthInMap(a);
0103     end
0104 end
0105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0106 
0107 %%%%%%%%% Cofactors not included in the analyses %%%%%%%%%%%%%%%%%%%%%%%%
0108 %%%%%%%%% Compartments included in the analyses %%%%%%%%%%%%%%%%%%%%%%%%%
0109 
0110 % propossed cofactors list for PART A
0111 cofactorsListA ={'h'; 'h2o'; 'nad'; 'nadh'; 'nadp'; 'nadph';...
0112     'fad'; 'fadh'; 'ppi'; 'adp'; 'gtp'; 'gdp'; 'gmp'; ...
0113     'utp'; 'udp'; 'ump'; 'ctp'; 'cdp'; 'cmp'};
0114 
0115 % propossed cofactors list for PART B
0116 cofactorsListB ={'h'; 'h2o'; 'co2'; 'coa'; 'accoa'; 'nad'; 'nadh'; 'nadp'; 'nadph';...
0117     'fad'; 'fadh'; 'pi'; 'ppi'; 'atp'; 'adp'; 'amp'; 'gtp'; 'gdp'; 'gmp';...
0118     'utp'; 'udp'; 'ump'; 'ctp'; 'cdp'; 'cmp'; 'gthox'; 'gthrd_c'; 'amet_c';...
0119     'ahcys_c'; 'glu__L'; 'gln__L'; 'akg'; 'nh4'};
0120 
0121 % compartments list
0122 compartmentsList = {'_e'; '_c'; '_m'; '_n'; '_x'; '_v'; '_r'; '_g'};
0123 
0124 % include compartment to cofactor lists A and B
0125 cofactorsListWcompA = cell(length(cofactorsListA), 1);
0126 c = 1; % general counter
0127 for i = 1:length(cofactorsListA)
0128     for j = 1:length(compartmentsList)
0129         cofactorsListWcompA{c, 1} = horzcat(cofactorsListA{i},compartmentsList{j});
0130         c = c + 1;
0131     end
0132 end
0133 cofactorsListWcompB = cell(length(cofactorsListB), 1);
0134 c = 1;
0135 for i = 1:length(cofactorsListB)
0136     for j = 1:length(compartmentsList)
0137         cofactorsListWcompB{c, 1}= horzcat(cofactorsListB{i},compartmentsList{j});
0138         c = c + 1;
0139     end
0140 end
0141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0142 
0143 
0144 %%%%%%%%%%%%%%%%%%%%%%%%%%%% PART A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145 
0146 % getting transport reactions
0147 [nrxns,~] = size(model2.rxns);
0148 transRxns = cell(nrxns, 1);
0149 for j = 1:size(model2.rxns)
0150     if true(contains(model2.subSystems {j}, 'S_Transport'))
0151         transRxns{j, 1} = model2.rxns{j};
0152     end
0153 end
0154 transRxns = transRxns(~cellfun('isempty',transRxns));
0155 
0156 % getting subStystems without transport
0157 SubSystems = (unique(getModelSubSystems(model2)));
0158 [nSubSystems, ~] = size(SubSystems);
0159 subSystemsNoTrans = cell(nSubSystems, 1);
0160 for i = 1:nSubSystems
0161     if ~contains(SubSystems{i,1},'S_Transport')&& ...
0162             ~contains(SubSystems{i, 1},'Extracellular exchange') &&...
0163             ~contains(SubSystems{i, 1},'S_tRNA_charging')
0164         subSystemsNoTrans{i, 1} = SubSystems{i, 1};
0165     else
0166         subSystemsNoTrans{i, 1} = {};
0167     end
0168 end
0169 subSystemsNoTrans = subSystemsNoTrans(~cellfun('isempty',subSystemsNoTrans));
0170 
0171 % defining subSystems of trasnport reactions to be colored, these
0172 % transport reactions complete reaction conections in subSystems which
0173 % require conexions among diferents compartments
0174 [nTransRxns,~] = size(transRxns);
0175 [nSubSystemsNoTrans, ~] = size (subSystemsNoTrans);
0176 totalRxnsToIncludeA = cell(nTransRxns,1);
0177 for i = 1:nSubSystemsNoTrans
0178     subSystemRxns = findRxnsFromSubSystem(model2, subSystemsNoTrans(i));
0179     metsSubSystemRxns = findMetsFromRxns(model2, subSystemRxns);
0180     [nMetsSubSystemRxns, ~] = size (metsSubSystemRxns);
0181     % getting unique compartments
0182     metsCompartment = cell(nMetsSubSystemRxns,1);
0183     for j = 1:nMetsSubSystemRxns
0184         str = metsSubSystemRxns{j};
0185         [~, s] = size(str);
0186         metsCompartment{j,1} = str(s-1:s);
0187     end
0188     metsCompartment = (unique(metsCompartment));
0189     [nUniqueMetsComp, ~] = size(metsCompartment);
0190     
0191     if nUniqueMetsComp > 1
0192         % getting duplicated mets like 'pyr_c' and 'pyr_m' which are 'pyr'
0193         % and pyr after erasing the compartment signature
0194         metsFromSubSystemRxnsNoComp = cell(nMetsSubSystemRxns,1);
0195         for j = 1:nMetsSubSystemRxns
0196             str = metsSubSystemRxns{j};
0197             [z, s] = size(str);
0198             metsFromSubSystemRxnsNoComp{j, 1} = str(1:s-2); % erasing comp signature
0199         end
0200         [~, o, ~] = unique(metsFromSubSystemRxnsNoComp,'first');
0201         indexToDupes = find(not(ismember(1:numel(metsFromSubSystemRxnsNoComp),o)));
0202         [~, nIndexToDupes] = size(indexToDupes);
0203         duplicatedMets = cell(nIndexToDupes,1);
0204         for j = 1:nIndexToDupes
0205             duplicatedMets(j, 1) = metsFromSubSystemRxnsNoComp (indexToDupes(1, j));
0206         end
0207         
0208         % getting reactions from duplicated mets
0209         dulicatedRxnMets = cell(nIndexToDupes * nUniqueMetsComp,1);
0210         c = 1; % duplicatedRxnMets counter
0211         for j = 1:nIndexToDupes
0212             str = duplicatedMets{j, 1};
0213             if ~ismember(str, cofactorsListA)
0214                 for k = 1:nUniqueMetsComp
0215                     str1 = metsCompartment{k};
0216                     if ~ismember(str1, compartmentsList(2)) % compartment '_c'
0217                         dulicatedRxnMets{c, 1} = {horzcat(str, '_c'), horzcat(str, str1)};
0218                         c = c + 1;
0219                     end
0220                 end
0221             end
0222         end
0223         dulicatedRxnMets = dulicatedRxnMets(~cellfun('isempty',dulicatedRxnMets));
0224         [ndulicatedRxnMets, ~] = size(dulicatedRxnMets);
0225         % from duplicated reactions, get reactions to subscribe a reactions
0226         % subsystem from non transport subsystems
0227         rxnsToInclude = cell(ndulicatedRxnMets, 1);
0228         for j = 1:ndulicatedRxnMets
0229             metsToFindInTrasnportRxns = dulicatedRxnMets{j, 1};
0230             for k = 1:nTransRxns
0231                 rxnFormula = printRxnFormula(model2, transRxns{k, 1}, false);
0232                 if true(contains(string(rxnFormula), metsToFindInTrasnportRxns (1, 1)))
0233                     if true(contains(rxnFormula, metsToFindInTrasnportRxns (1, 2)))
0234                         rxnsToInclude{j, 1} = transRxns{k, 1};
0235                     end
0236                 end
0237             end
0238         end
0239         rxnsToInclude = rxnsToInclude(~cellfun('isempty',rxnsToInclude));
0240         rxnsToInclude = unique(rxnsToInclude);
0241         [s, ~] = size(rxnsToInclude);
0242         totalRxnsToIncludeA = [totalRxnsToIncludeA; rxnsToInclude];
0243         
0244         for j = 1:s
0245             rxnID = findRxnIDs(model2, rxnsToInclude{j});
0246             model2.subSystems(rxnID) = {subSystemsNoTrans(i)};
0247         end
0248         
0249     end
0250 end
0251 totalRxnsToIncludeA = totalRxnsToIncludeA(~cellfun('isempty',totalRxnsToIncludeA));
0252 totalRxnsToIncludeA = unique(totalRxnsToIncludeA);
0253 
0254 %%%%%%%%%%%%%%%%%%%%%%%%%%%% PART B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0255 
0256 % Get subsystems category. This is based on the number of conexions that a
0257 % subsystem includes in the genome scale metabolic model2
0258 subsLength ={};
0259 for i = 1:length(subSystemsNoTrans)
0260     rxnsSubs = findRxnsFromSubSystem(model2,subSystemsNoTrans(i));
0261     metsSubsRxns = findMetsFromRxns(model2, rxnsSubs);
0262     metsSubsRxnsNoCof = {};
0263     for j = 1:length(metsSubsRxns)
0264         if ~ismember(metsSubsRxns(j),cofactorsListWcompB)
0265             metsSubsRxnsNoCof{j,1} = char(metsSubsRxns(j));
0266         end
0267     end
0268     metsSubsRxnsNoCof = metsSubsRxnsNoCof(~cellfun('isempty',metsSubsRxnsNoCof));
0269     
0270     metsRxns = findRxnsFromMets(model2, metsSubsRxnsNoCof);
0271     metsRxnsSubs = {};
0272     for j = 1:length(metsRxns)
0273         rxnID = findRxnIDs(model2,metsRxns(j));
0274         metsRxnsSubs{j,1} = string(model2.subSystems(rxnID));
0275     end
0276     metsRxnsSubs = unique(string(metsRxnsSubs));
0277     subsLength{i,1} = length(metsRxnsSubs);
0278     subsLength{i,2} = subSystemsNoTrans{i};
0279 end
0280 
0281 % Get  remaining transport reacctions, namely, reactions not included in
0282 % PART A
0283 
0284 [nTransRxns, ~] = size(transRxns);
0285 remainingTransRxns = cell(nTransRxns, 1);
0286 for i = 1:nTransRxns
0287     if ~ismember(transRxns(i),totalRxnsToIncludeA)
0288         remainingTransRxns{i,1} = transRxns{i, 1};
0289     end
0290 end
0291 remainingTransRxns = remainingTransRxns(~cellfun('isempty',remainingTransRxns));
0292 
0293 % Assign remaining transport reactions
0294 nonIncludSubs = {'S_tRNA_charging', 'S_Other',...
0295     'Biomass and maintenance functions', };
0296 totalRxnsToIncludeAB = totalRxnsToIncludeA;
0297 remTransRxnsWnewSubs = {};
0298 remTransRxnsSubs = {};
0299 rxnsToInclude = cell(length(remainingTransRxns), 1);
0300 subSToInclude = cell(length(remainingTransRxns), 1);
0301 for i = 1:length(remainingTransRxns)
0302     metsRemTransRxns = findMetsFromRxns(model2, remainingTransRxns(i));
0303     metsRxns = {};
0304     for j = 1:length(metsRemTransRxns)
0305         if true(contains(metsRemTransRxns{j},'_c') &&...
0306                 ~ismember(metsRemTransRxns(j), cofactorsListWcompB))
0307             metsRxns = findRxnsFromMets(model2,metsRemTransRxns(j));
0308         end
0309     end
0310     metsRxns = metsRxns(~cellfun('isempty',metsRxns));
0311     
0312     rxnSubs = {};
0313     for j = 1:length(metsRxns)
0314         rxnID = findRxnIDs(model2, metsRxns(j));
0315         rxnSubs{j,1} = model2.subSystems{rxnID};
0316     end
0317     rxnSubs = rxnSubs(~cellfun('isempty',rxnSubs));
0318     
0319     r = [];
0320     rankMax = 0;
0321     subSToUseStr = cell(1,1);
0322     for j = 1:length(rxnSubs)
0323         if true(~ismember(string(metsRxns(j)),totalRxnsToIncludeA))
0324             for k = 1:length(subsLength)
0325                 if true(and(ismember(rxnSubs{j}, subsLength{k,2}),...
0326                         ~ismember(rxnSubs{j}, nonIncludSubs)))
0327                     r(k, 1) = subsLength{k, 1};
0328                     subSToUseStr{1, 1} = string(subsLength(k,2));
0329                     if k > 1
0330                         if true(rankMax < r(k, 1))
0331                             rankMax = r(k, 1);
0332                             subSToUseStr{1, 1} = string(subsLength(k,2));
0333                         end
0334                     end
0335                 end
0336             end
0337         end
0338     end
0339     if ~isempty(subSToUseStr(1,1))
0340         rxnsToInclude{i} = remainingTransRxns{i};
0341         subSToInclude{i} = char(subSToUseStr{1, 1});
0342         remTransRxnsWnewSubs{i,1} = remainingTransRxns(i);
0343         remTransRxnsSubs{i,1} = subSToUseStr(1, 1);
0344         totalRxnsToIncludeAB = [totalRxnsToIncludeAB; rxnsToInclude];
0345     end
0346 end
0347 totalRxnsToIncludeAB = totalRxnsToIncludeAB(~cellfun('isempty',totalRxnsToIncludeAB));
0348 totalRxnsToIncludeAB = unique(totalRxnsToIncludeAB);
0349 
0350 for j = 1:length(rxnsToInclude)
0351     rxnID = findRxnIDs(model2, rxnsToInclude{j});
0352     model2.subSystems(rxnID) = {subSToInclude(j)};
0353 end
0354 
0355 %%%%%%%%%%%%%%%%%%%%%%%% ASIGN SUBSYSTEMS' COLOR %%%%%%%%%%%%%%%%%%%%%%%%%%%
0356 
0357 % Selection of reactions to include in the map
0358 if rxnsFluxTask == 1
0359     for i = 1:length(model2.rxns)
0360         if absFlxs(i)==0 % zero flux reactions are not included
0361             model2.subSystems(i) = {{'NonIncluded'}};  % included at the end
0362         end
0363     end
0364 elseif rxnsFluxTask == 2
0365     for i = 1:length(model2.rxns)
0366         if absFlxs(i) > 0.95 * max(absFlxs)% reactions causing growth (i.e. obj function) defect when deleting
0367             model2.subSystems(i) = {{'NonIncluded'}};   % included at the end
0368         end
0369     end
0370 elseif rxnsFluxTask == 3
0371     for i = 1:length(model2.rxns)
0372         if absFlxs(i) > 0.0 % zero flux reactions only, i.e. essential reations
0373             model2.subSystems(i) = {{'NonIncluded'}};   % included at the end
0374         end
0375     end
0376 elseif rxnsFluxTask == 4
0377     for i = 1:length(model2.rxns)
0378         if absFlxs(i) < 0.95 * max(absFlxs)% reactions fluxes causing no growth (i.e. obj function) defect when deleting
0379             model2.subSystems(i) = {{'NonIncluded'}};   % included at the end
0380         end
0381     end
0382 elseif rxnsFluxTask == 5
0383     for i = 1:length(model2.rxns)
0384         if absFlxs(i) < cutValue % zero flux reactions are not included
0385             model2.subSystems(i) = {{'NonIncluded'}};  % included at the end
0386         end
0387     end
0388 end
0389 
0390 %See colors letter-code at:
0391 %la.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/24497/versions/2/screenshot.PNG
0392 %[xmlMap, mapMap] = transformXML2Map('YeastMapUnified1.0b.xml');
0393 centralPathsList = {'S_GlycolysisGluconeogenesis', 'S_Pentose_Phosphate_Pathway',...
0394     'S_Pyruvate_Metabolism', 'S_Citric_Acid_Cycle', 'S_Taurine_Metabolism',...
0395     'S_Starch_and_Sucrose_Metabolism', 'S_Xylose_Metabolism',...
0396     'S_Fructose_and_Mannose_Metabolism', 'S_Galactose_metabolism',...
0397     'S_Arabinose_Metabolism', 'S_Alternate_Carbon_Metabolism', 'S_Methane_Metabolism'};
0398 aminoacPathsList = {'S_Glutamate_metabolism', 'S_Glutamine_Metabolism',...
0399     'S_Alanine_and_Aspartate_Metabolism', 'S_Histidine_Metabolism',...
0400     'S_Threonine_and_Lysine_Metabolism', 'S_Glycine_and_Serine_Metabolism',...
0401     'S_Tyrosine__Tryptophan__and_Phenylalanine_Metabolism',...
0402     'S_Arginine_and_Proline_Metabolism', 'S_Cysteine_Metabolism',...
0403     'S_Methionine_Metabolism', 'S_Valine__Leucine__and_Isoleucine_Metabolism',...
0404     'S_Other_Amino_Acid_Metabolism', 'S_Asparagine_metabolism', 'S_Nitrogen_Metabolism'};
0405 cofactorsPathsList = {'S_NAD_Biosynthesis', 'S_Pantothenate_and_CoA_Biosynthesis',...
0406     'S_Folate_Metabolism', 'S_Riboflavin_Metabolism', 'S_Quinone_Biosynthesis',...
0407     'S_Porphyrin_and_Chlorophyll_Metabolism'};
0408 lipidsPathsList = {'S_Fatty_Acid__Biosynthesis', 'S_Fatty_Acid_Degradation',...
0409     'S_Fatty_Acid_Metabolism', 'S_Glycerolipid_Metabolism', 'S_Glycoprotein_Metabolism',...
0410     'S_Phospholipid_Biosynthesis', 'S_Phospholipid_Metabolism',...
0411     'S_Sphingolipid_Metabolism', 'S_Sterol_Metabolism'};
0412 nucleotidsPathsList = {'S_Purine_and_Pyrimidine_Biosynthesis',...
0413     'S_Nucleotide_Salvage_Pathway', 'S_Thiamine_Metabolism', 'S_Pyridoxine_Metabolism'};
0414 enrgyPathsList = {'Biomass and maintenance functions', 'S_Anaplerotic_reactions',...
0415     'S_Oxidative_Phosphorylation'};
0416 colors = createColorsMap;
0417 
0418 if isempty(rxnLineColor)||rxnLineColor == "Pathways"
0419     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_GlycolysisGluconeogenesis', 'DARKBLUE');
0420     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pentose_Phosphate_Pathway', 'STEELBLUE');
0421     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pyruvate_Metabolism', 'DODGERBLUE');
0422     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Citric_Acid_Cycle', 'MEDIUMBLUE');
0423     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Taurine_Metabolism', 'MEDIUMBLUE');
0424     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Starch_and_Sucrose_Metabolism', 'BLUE');
0425     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Xylose_Metabolism', 'STEELBLUE');
0426     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fructose_and_Mannose_Metabolism', 'DODGERBLUE');
0427     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Galactose_metabolism', 'MEDIUMBLUE');
0428     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Arabinose_Metabolism', 'BLUE');
0429     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Alternate_Carbon_Metabolism', 'DARKBLUE');
0430     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Methane_Metabolism', 'DARKBLUE');
0431     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Taurine_Metabolism', 'DARKBLUE');
0432     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Anaplerotic_reactions', 'CRIMSON');
0433     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Complex_Alcohol_Metabolism', 'MAROON');
0434     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Oxidative_Phosphorylation', 'RED');
0435     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Purine_and_Pyrimidine_Biosynthesis', 'PINK');
0436     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Nucleotide_Salvage_Pathway', 'HOTPINK');
0437     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Thiamine_Metabolism', 'MAGENTA');
0438     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pyridoxine_Metabolism', 'PINK');
0439     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_NAD_Biosynthesis', 'PLUM');
0440     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Pantothenate_and_CoA_Biosynthesis', 'VIOLET');
0441     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Folate_Metabolism', 'DARKORCHID');
0442     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Riboflavin_Metabolism', 'PURPLE');
0443     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Quinone_Biosynthesis', 'MEDIUMORCHID');
0444     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Porphyrin_and_Chlorophyll_Metabolism', 'PLUM');
0445     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glutamate_metabolism', 'ORANGE');
0446     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glutamine_Metabolism', 'DARKORANGE');
0447     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Alanine_and_Aspartate_Metabolism', 'CORAL');
0448     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Histidine_Metabolism', 'GOLD');
0449     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Threonine_and_Lysine_Metabolism', 'ORANGE');
0450     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glycine_and_Serine_Metabolism', 'DARKORANGE');
0451     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Tyrosine__Tryptophan__and_Phenylalanine_Metabolism', 'CORAL');
0452     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Arginine_and_Proline_Metabolism', 'GOLD');
0453     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Cysteine_Metabolism', 'ORANGE');
0454     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Methionine_Metabolism', 'DARKORANGE');
0455     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Valine__Leucine__and_Isoleucine_Metabolism', 'CORAL');
0456     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Other_Amino_Acid_Metabolism', 'GOLD');
0457     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Asparagine_metabolism', 'ORANGE');
0458     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Nitrogen_Metabolism', 'DARKORANGE');
0459     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fatty_Acid__Biosynthesis', 'FORESTGREEN');
0460     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fatty_Acid_Degradation', 'MEDIUMSEAGREEN');
0461     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Fatty_Acid_Metabolism', 'DARKGREEN');
0462     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glycerolipid_Metabolism', 'SEAGREEN');
0463     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Glycoprotein_Metabolism', 'FORESTGREEN');
0464     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Phospholipid_Biosynthesis', 'MEDIUMSEAGREEN');
0465     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Phospholipid_Metabolism', 'DARKGREEN');
0466     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Sphingolipid_Metabolism', 'SEAGREEN');
0467     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Sterol_Metabolism', 'FORESTGREEN');
0468     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Other', 'BLACK');
0469     modifiedMap = colorSubsystem(modifiedMap, model2, 'Biomass and maintenance functions', 'RED');
0470 elseif rxnLineColor == "Subsystems"
0471     modifiedMap = colorSubsystem(modifiedMap, model2, centralPathsList, 'BLUE'); %mapMap to fluxMap
0472     modifiedMap = colorSubsystem(modifiedMap, model2, aminoacPathsList, 'ORANGE');
0473     modifiedMap = colorSubsystem(modifiedMap, model2, cofactorsPathsList, 'ORCHID');
0474     modifiedMap = colorSubsystem(modifiedMap, model2, lipidsPathsList, 'LIME');
0475     modifiedMap = colorSubsystem(modifiedMap, model2, nucleotidsPathsList, 'MAGENTA');
0476     modifiedMap = colorSubsystem(modifiedMap, model2, enrgyPathsList, 'RED');
0477     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Complex_Alcohol_Metabolism', 'GOLDENROD');
0478     modifiedMap = colorSubsystem(modifiedMap, model2, 'S_Other', 'GOLDENROD');
0479 elseif ~(rxnLineColor == "Subsystems" || rxnLineColor == "Payhways" || isempty(rxnLineColor))
0480     for i = 1:length(modifiedMap.rxnName)
0481         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0482         if ~isempty(a) & string(model2.subSystems(a)) ~= "NonIncluded"
0483             modifiedMap.rxnColor(i, 1) = {colors(rxnLineColor)};
0484         end
0485     end
0486 end
0487 if rxnLineColor == "Subsystems" || rxnLineColor == "Payhways" || isempty(rxnLineColor)
0488     h2o = {'H2Ot', 'H2Otm',    'H2Otn', 'H2Otv', 'H2Otp', 'H2Oter'};
0489     for i = 1:length(modifiedMap.rxnName)
0490         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0491         for j = 1:length(h2o)
0492             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(h2o(j))...
0493                     & string(model2.subSystems(a)) ~= "NonIncluded"
0494                 modifiedMap.rxnColor(i, 1) = {colors('CYAN')};
0495             end
0496         end
0497     end
0498     o2 = {'O2t', 'O2tm', 'O2ter', 'O2tp'};
0499     for i = 1:length(modifiedMap.rxnName)
0500         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0501         for j = 1:length(o2)
0502             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(o2(j))...
0503                     & string(model2.subSystems(a)) ~= "NonIncluded"
0504                 modifiedMap.rxnColor(i, 1) = {colors('CYAN')};
0505             end
0506         end
0507     end
0508     Pi = {'PIt2r', 'PIt2m', 'PIt5m', 'PIt2n', 'PIt2v', 'PIt2p', 'ATPtm_H', 'ATPtp_H', ...
0509         'ATP2tp_H', 'ATPS3g', 'ATPS3v', 'UTPtm', 'CTPtm', 'GTPt2m', 'ATPS'};
0510     for i = 1:length(modifiedMap.rxnName)
0511         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0512         for j = 1:length(Pi)
0513             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(Pi(j))...
0514                     & string(model2.subSystems(a)) ~= "NonIncluded"
0515                 modifiedMap.rxnColor(i, 1) = {colors('RED')};
0516             end
0517         end
0518     end
0519     co2 = {'CO2t', 'CO2tg', 'CO2tm', 'CO2tn', 'HCO3tn', 'CO2tv', 'CO2tp'};
0520     for i = 1:length(modifiedMap.rxnName)
0521         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0522         for j = 1:length(co2)
0523             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(co2(j))...
0524                     & string(model2.subSystems(a)) ~= "NonIncluded"
0525                 modifiedMap.rxnColor(i, 1) = {colors('SLATEGRAY')};
0526             end
0527         end
0528     end
0529     aa = {'MMETt2', 'GLUt2r', 'GLNt2r', 'GLNt7', 'GLUt7', 'LCYSTintv', 'NH4t', 'TCHOLAabcv'};
0530     for i = 1:length(modifiedMap.rxnName)
0531         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0532         for j = 1:length(aa)
0533             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(aa(j))...
0534                     & string(model2.subSystems(a)) ~= "NonIncluded"
0535                 modifiedMap.rxnColor(i, 1) = {colors('ORANGE')};
0536             end
0537         end
0538     end
0539     cenMet = {'COAtim', 'NAt3_1', 'MALS'};
0540     for i = 1:length(modifiedMap.rxnName)
0541         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0542         for j = 1:length(cenMet)
0543             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(cenMet(j))...
0544                     & string(model2.subSystems(a)) ~= "NonIncluded"
0545                 modifiedMap.rxnColor(i, 1) = {colors('BLUE')};
0546             end
0547         end
0548     end
0549     lip = {'G3PIt'};
0550     for i = 1:length(modifiedMap.rxnName)
0551         a = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0552         for j = 1:length(lip)
0553             if ~isempty(a) & string(modifiedMap.rxnName{i}) == string(lip(j))...
0554                     & string(model2.subSystems(a)) ~= "NonIncluded"
0555                 modifiedMap.rxnColor(i, 1) = {colors('LIME')};
0556             end
0557         end
0558     end
0559 end
0560 
0561 if rxnsFluxTask == 5 % change color to dark-color in negative values
0562     for i = 1:length(modifiedMap.rxnName)
0563         valueskeys = values(colors)';
0564         valueskeys (:,2) = keys(colors)';
0565         a = find(ismember(valueskeys, modifiedMap.rxnColor{i, 1}));
0566         b = find(ismember(model2.rxns, modifiedMap.rxnName{i}));
0567         if ~isempty(b) & v1(b) < 0 & string(model2.subSystems(b)) ~= "NonIncluded"
0568             if string(valueskeys (a,2)) == 'LIME'
0569                 modifiedMap.rxnColor(i, 1) = {colors('GREEN')};
0570             else
0571                 rxnColor =strcat('DARK', char(valueskeys (a(1,1),2)));
0572                 if string(rxnColor) == 'DARKAQUA'
0573                     rxnColor = 'DARKCYAN';
0574                 elseif string(rxnColor) == 'DARKFUCHSIA'
0575                     rxnColor = 'DARKMAGENTA';
0576                 end
0577                 modifiedMap.rxnColor(i, 1) = {colors(rxnColor)};
0578             end
0579         end
0580     end
0581 end
0582 
0583 end