ecFVA

 ecFVA
   Flux variability analysis is performed on the ecModel, and isozymic
   reactions are combined to construct ouput minFlux and maxFlux vectors,
   which follow the same order of model.rxns. The output from this
   function does not include enzyme usage reactions, to observe these, on
   could consider running flux variability directly on the ecModel.

 Input:
   ecModel     an ecModel in GECKO 3 format (with ecModel.ec structure)
   model       non-ecModel variant of the ecModel, to which the minFlux
               and maxFlux will be mapped
 Output:
   minFlux     vector of minimum flux rates, corresponding to model.rxns
   maxFlux     vector of maximum flux rates, corresponding to model.rxns

 Usage: [minFlux, maxFlux] = ecFVA(ecModel, model)

0001 function [minFlux, maxFlux] = ecFVA(ecModel, model)
0002 % ecFVA
0003 %   Flux variability analysis is performed on the ecModel, and isozymic
0004 %   reactions are combined to construct ouput minFlux and maxFlux vectors,
0005 %   which follow the same order of model.rxns. The output from this
0006 %   function does not include enzyme usage reactions, to observe these, on
0007 %   could consider running flux variability directly on the ecModel.
0008 %
0009 % Input:
0010 %   ecModel     an ecModel in GECKO 3 format (with ecModel.ec structure)
0011 %   model       non-ecModel variant of the ecModel, to which the minFlux
0012 %               and maxFlux will be mapped
0013 % Output:
0014 %   minFlux     vector of minimum flux rates, corresponding to model.rxns
0015 %   maxFlux     vector of maximum flux rates, corresponding to model.rxns
0016 %
0017 % Usage: [minFlux, maxFlux] = ecFVA(ecModel, model)
0018 
0019 rxnIDs = regexprep(ecModel.rxns,'(_REV)?(_EXP_\d+)?','');
0020 [rxnIDmap, convRxnID] = findgroups(rxnIDs);
0021 
0022 solMaxAll = nan(numel(ecModel.rxns),numel(convRxnID));
0023 solMinAll = solMaxAll;
0024 
0025 pool = gcp('nocreate');
0026 if isempty(pool)
0027     parpool;
0028 end
0029 
0030 D = parallel.pool.DataQueue;
0031 progressbar('Running ecFVA');
0032 afterEach(D, @nUpdateProgressbar);
0033 
0034 N = numel(convRxnID);
0035 p = 1;
0036 
0037 parfor i=1:N
0038     tmpModel = ecModel;
0039     tmpModel.c = zeros(numel(tmpModel.rxns),1);
0040 
0041     rxnsToOptim = find(rxnIDmap == i);
0042     rxnsOptIDs  = ecModel.rxns(rxnsToOptim);
0043     rxnsToMin = endsWith(rxnsOptIDs,'_REV') | contains(rxnsOptIDs,'_REV_EXP_');
0044     rxnsToMax = rxnsToOptim(~rxnsToMin);
0045     rxnsToMin = rxnsToOptim(rxnsToMin);
0046 
0047     tmpModel.c(rxnsToMax) = 1;
0048     tmpModel.c(rxnsToMin) = -1;
0049     solMax=solveLP(tmpModel);
0050     if ~isempty(solMax.x)
0051         solMaxAll(:,i)=solMax.x;
0052     end
0053     tmpModel.c(rxnsToMax) = -1;
0054     tmpModel.c(rxnsToMin) = 1;
0055     solMin=solveLP(tmpModel);
0056     if ~isempty(solMin.x)
0057         solMinAll(:,i)=solMin.x;
0058     end    
0059     send(D, i);
0060 end
0061 
0062 minFlux=min(solMinAll,[],2,'omitnan');
0063 maxFlux=max(solMaxAll,[],2,'omitnan');
0064 
0065 mappedFlux = mapRxnsToConv(ecModel,model,[minFlux maxFlux]);
0066 
0067 minFlux=mappedFlux(:,1);
0068 maxFlux=mappedFlux(:,2);
0069 
0070 % Mapped flux might have swapped directionality: min/max might be swapped
0071 swapDir = minFlux > maxFlux;
0072 if any(swapDir)
0073     tmpFlux = minFlux(swapDir);
0074     minFlux(swapDir) = maxFlux(swapDir);
0075     maxFlux(swapDir) = tmpFlux;
0076 end
0077 
0078 function nUpdateProgressbar(~)
0079 progressbar(p/N);
0080 p = p + 1;
0081 end
0082 end