ObjComparison/maxATPProd.m at master · rs-costa/ObjComparison · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
% Exploring the cellular objective in flux balance constraint-based models.

% Author: Rafael Costa
% Affiliation: Instituto de Engenharia de Sistemas e Computadores, Investigacão e Desenvolvimento (INESC-ID), Lisboa
% and
% Center for Intelligent Systems, LAETA, IDMEC, IST, University of Lisbon
% Author: Nguyen Hoang Son

function FBAsolution=maxATPProd(glob,model,netcode,taskcode,minNorm)
%Max ATP production
    if exist('minNorm','var')
        if isempty(minNorm)
            minNorm=true;
        end
    else
        minNorm=true;
    end

    if netcode==1
        [modelIrrev, matchRev, rev2irrev, irrev2rev]= convertToIrreversible(model);
        atpId=findMetIDs(modelIrrev,'atp[c]');
        bmId=findRxnIDs(modelIrrev,'Biomass_Ecoli_core_w_GAM');
    %     obj_rnxs={'PGK_f','PYK','SUCOAS_f','ATPS4r_f','ACKr_b'};

    elseif netcode==2
        [modelIrrev, matchRev, rev2irrev, irrev2rev]= convertToIrreversible(model);
        atpId=findMetIDs(modelIrrev,'ATP');
        bmId=findRxnIDs(modelIrrev,'biomass');
    %     obj_rnxs={'pgk_f','pykA','pykF','sucCD_f','atp'...
    %                'ackA_f','ackB_f','tdcD_f','purT_f'};
    elseif netcode==3
        [modelIrrev, matchRev, rev2irrev, irrev2rev]= convertToIrreversible(model);
        atpId=findMetIDs(modelIrrev,'atp[c]');
        bmId=findRxnIDs(modelIrrev,'Ec_biomass_iAF1260_core_59p81M');
    %     obj_rnxs={'PGK_f','PYK','SUCOAS_f','ATPS4rpp_f','ACKr_b'};
    end
    %%Get list of reactions which produce ATP
    obj_rnxs=full(modelIrrev.S);
    obj_rnxs=obj_rnxs(atpId,:);
    obj_rnxs(obj_rnxs<0)=0;
    obj_rnxs(bmId)=0;

    modelIrrev.csense=model.csense;
    modelIrrev.c=obj_rnxs';
    if minNorm
        FBAsolution = optimizeCbModel(modelIrrev,'max','one');
    else
        FBAsolution = optimizeCbModel(modelIrrev,'max');
    end

    if FBAsolution.stat < 1
        FBAsolution.x=[];
        FBAsolution.f=0;
        return;
    end
    x_irrev=zeros(length(rev2irrev),1);
    for i=1:length(rev2irrev)
        x_irrev(i)=sum(FBAsolution.x(rev2irrev{i}));
    end
    FBAsolution.x=x_irrev;
    if minNorm
        FBAsolution.minE=eclDistance(glob,FBAsolution.x)/100;
        FBAsolution.maxE=FBAsolution.minE;
    else
        A=[model.S(model.csense=='L',:);-model.S(model.csense=='G',:)];
        b=[model.b(model.csense=='L',:);-model.b(model.csense=='G',:)];
        Aeq=model.S(model.csense=='E',:);
        beq=model.b(model.csense=='E',:);

        x0=FBAsolution.x;
        f0=FBAsolution.f;

        Aeq=[Aeq;model.c'];
        beq=[beq;f0];

        [FBAsolution.minE,FBAsolution.maxE,FBAsolution.x]=fidelity4linear(glob,taskcode,x0,Aeq,beq,A,b,model.lb,model.ub);
    end

end