function [futileCycleReactions,FutileCyclesTest]=identifyFutileCycles(model)
% This function serves to debug a refined reconstruction generated by the
% DEMETER pipeline that produced infeasible amounts of ATP. The output 
% futileCycleReactions shows reactions that are involved in futile cycles.
% To eliminate the futile cycle, one of these reactions needs to be
% replaced by an irreversible version of the same reaction. To accomplish
% this, the replacement step needs to be added to the function 
% removeFutileCycles. It is requested that identified solutions are
% submitted as a pull request.
%
% USAGE:
%
%   [gapfilledReactions,replacedReactions,revisedModel]=debugModel(model,testResultsFolder, inputDataFolder,reconVersion,microbeID,biomassReaction)
%
% INPUT
% model:                   COBRA model structure
%
% OUTPUT
% futileCycleReactions:    Reactions that are involved in futile cycles and
%                          need to be constrained
% FutileCyclesTest         Computed fluxes for ATP production
%
% .. Authors:
%       - Almut Heinken, 2016-2020
    
% load complex medium
constraints = readtable('ComplexMedium.txt', 'Delimiter', 'tab');
constraints=table2cell(constraints);
constraints=cellstr(string(constraints));

% apply complex medium
model = useDiet(model,constraints);
model=changeRxnBounds(model,'EX_o2(e)',-10,'l');
model=changeObjective(model,'DM_atp_c_');
FBAorg=optimizeCbModel(model,'max');

if FBAorg.f > 150

Rxns=printRxnFormula(model);
for i=1:length(model.rxns)
    if ~strcmp(model.rxns{i},'DM_atp_c_') && ~strncmp(model.rxns{i},'ATPS',4)
        modelTest=model;
        FutileCyclesTest{i,1}=model.rxns{i,1};
        modelTest.lb(i,1)=0;
        modelTest.ub(i,1)=0;
        FBA=optimizeCbModel(modelTest,'max');
        FutileCyclesTest{i,2}=Rxns{i,1};
        FutileCyclesTest{i,3}=model.lb(i,1);
        FutileCyclesTest{i,4}=model.ub(i,1);
        FutileCyclesTest{i,5}=FBA.f;
        FutileCyclesTest{i,6}=FBAorg.x(i,1);
    end
end

FutileCyclesTest(find(cellfun(@isempty, FutileCyclesTest(:,1))),:)=[];
futileCycleReactions=FutileCyclesTest(cell2mat(FutileCyclesTest(:,5))<150,1);

else
    FutileCyclesTest = {};
    futileCycleReactions = {};
end

end
