function [Leaks, Leaks_min, Leaks_max, varargout]=bktr_LeakFVA(irrevGSM,GSM,matchRev,map,scale,rev2irrev, perc_opt)

if nargin==6
    perc_opt = 1;
end

%% Run normal FBA
MinimizedFlux = optimizeCbModel(irrevGSM,'min');

%% Prepare FVA
% Fix opt flux 
allowed_diff =abs( MinimizedFlux.f -(MinimizedFlux.f * perc_opt));
irrevGSM.lb(find(irrevGSM.c))=MinimizedFlux.f -allowed_diff; % will never be reached
irrevGSM.ub(find(irrevGSM.c))=MinimizedFlux.f +allowed_diff;

% select reactions involved in leaks (to save time)
mapped_mets=[map.GSMMets{:}];
rxns_of_int=findRxnsFromMets(GSM, GSM.mets(mapped_mets));
[a,b]= ismember(rxns_of_int,GSM.rxns);
rxns_of_int=rev2irrev(b(a));
exrn  = rev2irrev(find(findExcRxns(GSM)));

selected=unique([ [exrn{:}] [rxns_of_int{:}]]); %irrev model
rxns_of_int_orig =unique([ b(a)' find(findExcRxns(GSM))']); % reversible model

 
%% FVA
% initialize storers
vIrrev_min= NaN(size(irrevGSM.lb)); vIrrev_max= NaN(size(irrevGSM.ub));

fprintf('%4s\t%4s\t%10s\t%9s\t%9s\t%9s\n','No','Perc','Name','Min','OPT','Max');count=0;
for k=selected %run FVA on all selected reactions
    count=count+1;
    irrevGSM.c=zeros(size(irrevGSM.c));
    irrevGSM.c(k)=1;
    
    min_sol=optimizeCbModel(irrevGSM,'min');
    vIrrev_min(k)=min_sol.f;
    
    irrevGSM.c(k)=-1;
    max_sol=optimizeCbModel(irrevGSM,'min');
    vIrrev_max(k)=-max_sol.f;
    
    fprintf('%4d\t%4.0f\t%10s\t%9.4f\t%9.4f\t%9.4f\n',count,100*count/length(selected),...
        irrevGSM.rxns{k}, vIrrev_min(k),MinimizedFlux.x(k),vIrrev_max(k));
    
end

% Convert results to reversible form
vRev = convertIrrevFluxDistribution(MinimizedFlux.x,matchRev);
[vRev_min, vRev_max] = bktr_convertIrrevFluxDistributionFVABounds(vIrrev_min,vIrrev_max,matchRev);


%% Collect Leaks
% find Excretion rxns and their excr metabolites
IdxExcrRxns=find(findExcRxns(GSM,1));
[IdxExcrMets,~]=find(GSM.S(:,IdxExcrRxns));

   
% all mapped rxns
rxns_of_int=~cellfun(@ischar, map.GSMRxns);
rxns_of_int=unique([map.GSMRxns{rxns_of_int}]);
% add lumped rxns if applicable
if isfield(map,'GSMLumpedRxns')
    rxns_of_int = unique([rxns_of_int map.GSMLumpedRxns]);
end

[a,b]=ismember(rxns_of_int_orig,rxns_of_int);

for k=1:length(map.DMets)
    if ~ischar(map.GSMMets{k})
        % find Rxn excreting the metabolites mapped to the current dyn metabolite
        mapped_mets=[map.GSMMets{k}];
        
        % flux
        Sw=full(GSM.S).*repmat(vRev,1,size(GSM.S,1))';
        Sw(:,rxns_of_int)=[];
        Leaks(k)=sum(Sw(mapped_mets,:),2)*scale;
        
        % min flux
        Sw_min=full(GSM.S).*repmat(vRev_min,1,size(GSM.S,1))';
        Sw_min=Sw_min(:,rxns_of_int_orig(~a));
        
        % min flux
        Sw_max=full(GSM.S).*repmat(vRev_max,1,size(GSM.S,1))';
        Sw_max=Sw_max(:,rxns_of_int_orig(~a));
        
        % account for influx-outflux combinatorics
        for j=1: length(mapped_mets)
            Sw=full(GSM.S);
            Sw=Sw(:,rxns_of_int_orig(~a));
            
            max_vals=(Sw_max(mapped_mets(j),:));
            min_vals=(Sw_min(mapped_mets(j),:));
            
            signs = sign(Sw(mapped_mets(j),:));
            cons = signs<0;
            prod = signs>0;
            
            maxval(j) = sum(min_vals(cons)) + sum(max_vals(prod));
            minval(j) = sum(max_vals(cons)) + sum(min_vals(prod));
            
        end
        
        
        Leaks_max(k)=sum(maxval)*scale;
        Leaks_min(k)=sum(minval)*scale;
        
        
    else
        Leaks(k)=0;
    end
    
end
% Other required Uptakes (optional)
idx_others =  abs(vRev(IdxExcrRxns))>1e-10; % exclude numerical errors
AddUptakes =  table(GSM.rxnNames(IdxExcrRxns(idx_others)),vRev(IdxExcrRxns(idx_others))*scale,...
    vRev_min(idx_others)*scale, vRev_max(idx_others)*scale);


varargout={AddUptakes};
