sparseFBA Not Minimizing Number of Reactions

[Matthew Scarborough]

I am trying to use sparseFBA with a model I built with cobraPy. I upload the model to Matlab as an SBML using 

readCbModel()

I also verify the model with 

verifyModel()

I get one error because I do not have any GPRs or genes in the model, but the model seems to work for downstream analyses (except sparseFBA).

When I use sparseFBA, I get the same reaction fluxes I get with pFBA (in cobraPy). However, I can knock out a set of genes that reduces the number of reactions (to 34 from 42) while maintaining the objective function at its maximum (in my case, ATP hydrolysis, to within 13 decimal places). For sparseFBA I use: 

[vSparse, sparseRxnBool, essentialRxnBool] = sparseFBA(model,'max',1,1,'all',1)

A list of reaction fluxes for the "sparseFBA" solution and the sparseFBA solutions + knockouts are located in this google sheet.

Lastly, I have used both glpk and gurobi solvers which give the same solutions. I am using Matlab2016b.

I would expect sparseFBA to minimizing the total number of reactions with flux while maintaining the maximum objective function. Is this correct? If so, do you have any suggestions to help me get sparseFBA working correctly?

Thanks for the help, and please let me know if you need any additional information.

Matt

[Matthew Scarborough]

Here is the COBRAconfigReport

[Matthew Scarborough]

I also tried the method shown in the Sparse FBA Tutorial (https://opencobra.github.io/cobratoolbox/latest/tutorials/tutorialSparseFBA.html), and I got the following output and error:

>> feasTol = getCobraSolverParams('LP', 'feasTol')

feasTol =

   1.0000e-09

>> minInf = -1000;

maxInf = 1000;

printConstraints(model, minInf, maxInf);

MinConstraints:

EX_xyl__D_e -1

EX_ac_e -3

maxConstraints:

>> osenseStr='max';

>> minNorm='zero'

minNorm =

zero

>> sparseFBAsolution = optimizeCbModel(model, osenseStr, minNorm)

sparseFBAsolution = 

  struct with fields:

         full: [115×1 double]

          obj: 3.3333

        rcost: []

         dual: []

        slack: [92×1 double]

       solver: 'gurobi'

    algorithm: 'default'

         stat: 1

     origStat: 'OPTIMAL'

         time: 0.0400

        basis: [1×1 struct]

            f: 3.3333

            x: [115×1 double]

            v: [115×1 double]

            w: []

            y: []

            s: [92×1 double]

>> v = sparseFBAsolution.v

v =

         0

         0

         0

         0

         0

         0

   -1.0000

    1.0000

         0

    1.0000

    1.0000

         0

   -0.3333

   -0.3333

    0.3333

    0.3333

    0.3333

         0

         0

         0

         0

         0

         0

    0.6667

    0.6667

    0.6667

    1.6667

   -1.6667

   -1.6667

    1.6667

    1.6667

         0

         0

         0

         0

    0.5556

   -0.5556

         0

    1.6667

    0.5556

    0.5556

    0.5556

    0.5556

         0

         0

         0

         0

         0

    0.5556

    0.5556

    0.5556

    0.5556

         0

         0

         0

         0

         0

    0.5556

    0.5556

    0.5556

    0.5556

         0

         0

    0.5556

   -0.5556

    0.5556

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

    3.3333

    1.6667

   -1.6667

    1.6667

         0

         0

   -1.6667

    1.6667

    3.3333

>> nonZeroFlag = 1;

printFluxVector(model, v, nonZeroFlag);

EX_xyl__D_e                  -1

XYLt                           1

XYLI1                         1

XYLK                           1

RPE                     -0.3333

RPI                     -0.3333

TKT1                      0.3333

TALA                      0.3333

TKT2                      0.3333

PFK                      0.6667

FBA                      0.6667

TPI                      0.6667

GAPD                       1.667

PGK                      -1.667

PGM                      -1.667

ENO                       1.667

PYK                       1.667

ACtr                      0.5556

EX_ac_e                 -0.5556

PFOR                       1.667

ACACT1r                  0.5556

HACD1                    0.5556

ECOAH1                    0.5556

EBACD1                    0.5556

ACACT2                    0.5556

HACD2                    0.5556

ECOAH2                    0.5556

EBACD2                    0.5556

ACACT3                    0.5556

HACD3                    0.5556

ECOAH3                    0.5556

EBACD3                    0.5556

CoATC8                    0.5556

octat                   -0.5556

EX_octa_e                0.5556

RNF1                       3.333

ATPS4r                     1.667

H2Ot                      -1.667

EX_h2o_e                   1.667

co2t                      -1.667

EX_co2_e                   1.667

ATP_Hydrolysis             3.333

>> fprintf('%u%s\n',nnz(v),' active reactions in the sparse flux balance analysis solution.');

42 active reactions in the sparse flux balance analysis solution.

>> [c,S,b,lb,ub,csense] = deal(model.c,model.S,model.b,model.lb,model.ub,model.csense);

[m,n] = size(S);

LPproblem = struct('c',c,'osense',-1,'A',S,'csense',csense,'b',b,'lb',lb,'ub',ub);

>> LPsolution = solveCobraLP(LPproblem);

if LPsolution.stat == 1

vFBA = LPsolution.full(1:n);

else

vFBA = [];

error('FBA problem error!')

end

>> fprintf('%u%s\n',nnz(vFBA),' active reactions in the flux balance analysis solution.')

38 active reactions in the flux balance analysis solution.

>> approximations = {'cappedL1','exp','log','SCAD','lp-','lp+'};

>> constraint.A = [S ; c'];

constraint.b = [b ; c'*vFBA];

constraint.csense = [csense;'E'];

constraint.lb = lb;

constraint.ub = ub;

>> bestResult = n;

bestAprox = '';

for i=1:length(approximations)

solution = sparseLP(char(approximations(i)),constraint);

if solution.stat == 1

nnzSol=nnz(abs(solution.x)>feasTol);

fprintf('%u%s%s',nnzSol,' active reactions in the sparseFBA solution with ', char(approximations(i)));

if bestResult > nnzSol

bestResult=nnzSol;

bestAprox = char(approximations(i));

solutionL0 = solution;

end

end

end

Error using sparseLP (line 95)

Error:LHS matrix is not defined

On Friday, December 28, 2018 at 10:17:48 PM UTC-5, Matthew Scarborough wrote:

[Matthew Scarborough]

One more follow-up... when I run the tuorial_sparseFBA.m, I get the same error as above. 

Error using sparseLP (line 95)

Error:LHS matrix is not defined

This is all making me think there is some sort of compatibility issue with my version of Matlab and the solver. 

[Matthew Scarborough]

Sorry for all the postings, but here is the latest--

I reverted to an older version of MacOSX (10.13) and installed Matlab2016b with Gurobi 8.0.1. (a winning combination according to the compatibility table). The sparseFBA command is still not giving me the sparsest solution-- the solutions still match the fluxes here.

The updated config report is attached. 

Thanks in advance for your help with this issue! I look forward to your reply. 

[Ronan M.T. Fleming]

Hi Matt,

please send me your .mat file with the model in it. I can test it next week and see what is happening.

Regards,

Ronan

--

---
You received this message because you are subscribed to the Google Groups "COBRA Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbo...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

--
Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.
----------------------------------------------------------------------------
Assistant Professor,

Division of Systems Biomedicine and Pharmacology,

Leiden Academic Centre for Drug Research,

Faculty of Science,

Leiden University.

https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming

&

H2020 Project Coordinator

Systems Medicine of Mitochondrial Parkinson’s Disease

http://sysmedpd.eu

----------------------------------------------------------------------------
Mobile:  +353 873 413 072
Skype: ronan.fleming
----------------------------------------------------------------------------
(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)