about .summary() usage

[gcxxx]

Hi everyone:

I load the iCHOv1(from BIGG) model into cobrapy and use model.optimize() to get fba result. Then I use model.summary() to see the profile of metabolism , # subsequently, when I try to use model.reaction.flux , an error happen

Why could that happen?

[Midnighter]

Hi @gcxxx,

Instead of describing your steps in words, can you please provide the exact code that you used? Thank you.

[Midnighter]

Glad to be of help! It's unexpected that it doesn't work with CPLEX, though. Can you share your model so we can try to reproduce that?

[gcxxx]

@Midnighter
You can download them from this link : http://bigg.ucsd.edu/models/iCHOv1

[gcxxx]

@Midnighter Sorry to bother you again🤦‍，but it seems like .summary() have a bug?
cmodel = read_sbml_model('d:/data/chomodel/iCHOv1_final.xml') csol = cmodel.optimize() cmodel.metabolites.g3p_c.summary() cmodel.reactions.LDH_Lm.flux,csol.fluxes.loc['LDH_Lm']

cmodel = read_sbml_model('d:/data/chomodel/iCHOv1_final.xml') csol = cmodel.optimize() cmodel.reactions.LDH_Lm.flux,csol.fluxes.loc['LDH_Lm']
Two code block above should have the same reslt? But the result is abnormal

[gcxxx]

I use the glpk solver

[gcxxx]

It seems like .summary() usage will change the final result

[cdiener]

Not really a bug, but if you don't pass a solution to summary it will generate a new one and this may be different from a previous one. Note that FBA solutions are generally not unique. See #1029 for some workarounds.

[gcxxx]

@cdiener That is helpful. Thank you!

[gcxxx]

A little confusion about not unique solution:

I use the textbook model to test whether the fba solution is unique, I solve a hundred time but each time gets the same result. BTW, I use glpk solver.

[cdiener]

Do you mean the objective value? That one is guaranteed to be unique, but the individual fluxes in the model are not. You can use a flux variability analysis to see how much they may vary.

[gcxxx]

First of all,Thanks your help!

a = cobra.test.create_test_model('textbook')
a.optimize()

code above always give the same result(both the objective value and fluxes) I run the code a dozens of time.
1

Btw, I actually can not understand FVA exactly. Does FVA result can show me the fluxes range of each reaction when the objective value is maximum? What does fraction_of_optimum exactly mean? I have read the documents but I am still puzzled.It seems like 0.9 fraction_of_optimum will greatly impact the range.

[bdelepine]

Hi everyone!

I've run some tests, out of curiosity. I could reproduce the example of @gcxxx and I may understand the confusing part.

I downloaded iCHOv1 from http://bigg.ucsd.edu/models/iCHOv1.

First, let's see what the code from this reply gives us:

cmodel = cobra.io.read_sbml_model("iCHOv1.xml")
cmodel.solver = "cplex"  # case 1
cmodel.solver = "glpk"   # case 2

csol = cmodel.optimize()
print(cmodel.metabolites.g3p_c.summary())  # we use .summary()
csol.status

In both cases, we get something along:

g3p_c
=====
Formula: C3H5O6P

Producing Reactions
-------------------
Percent    Flux Reaction                          Definition
 37.48% 0.06743      FBA            fdp_c <=> dhap_c + g3p_c
  9.85% 0.01772     TKT1 r5p_c + xu5p__D_c <=> g3p_c + s7p_c
  9.85% 0.01772     TKT2 e4p_c + xu5p__D_c <=> f6p_c + g3p_c
 42.83% 0.07706      TPI                    dhap_c <=> g3p_c

Consuming Reactions
-------------------
Percent    Flux Reaction                                      Definition
100.00% -0.1799     GAPD g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c

'optimal'

However the result of the line cmodel.reactions.LDH_Lm.flux is different if we use CPLEX or GLPK:

• CPLEX yields a CplexSolverError and at the root of the traceback we have OptimizationError: Likely no solution exists. Original solver message: CPLEX Error 1217: No solution exists...
• GLPK simply get back a flux of 0.0.

So what happened? My guess is that LDH_Lm cannot carry any flux and that GLPX rounds it to 0 while CPLEX explicitly warns that it cannot carry flux. Let me note here that accessing the flux value from the reaction object is considered bad practice as stated in the documentation. You should always access fluxes from a Solution object... and it where the fun begins!

The line csol.fluxes.loc['LDH_Lm'] returns in both cases 0.0. My guess is that cobra is smart enough to catch the (eventual) error and displays the expected value of 0.0... but I understand it is confusing since accessing "directly" the flux from the reaction object (which is bad!) may yield an error.

And then comes the .summary() part for even more fun:

cmodel = cobra.io.read_sbml_model("iCHOv1.xml") 
cmodel.solver = "cplex"  # case 1
cmodel.solver = "glpk"   # case 2

# first,  we use .summary()
csol = cmodel.optimize()
print(cmodel.metabolites.g3p_c.summary())

# second, we do not use .summary()
csol = cmodel.optimize()
cmodel.reactions.LDH_Lm.flux, csol.fluxes.loc['LDH_Lm']

I get:

• CPLEX: (0.16516256421589862, 0.16516256421589862)
• GLPK: (5.551115123125783e-17, 5.551115123125783e-17)
• commenting out the .summary() line, I get (0.0, 0.0) with both solvers.

This is certainly confusing... but not unexpected. So what is happening? First, let's note that 1e-17 is 0 considering the solver's precision. But what about the 0.16 flux value given by CPLEX? This is where you need to understand what @cdiener said: only the objective value is guaranteed to be optimal, but several flux distributions may yield the same optimal objective value. This is fairly common (and the sign that the model is not constrained enough). Only one flux distribution is reported when you do an FBA. If you wish to get the minimal and maximal bound of each of your reactions, well, you can ask the solver for precisely that:

1. first get the optimal FBA objective value
2. then, for each reaction, constrain the model's previous objective to its optimal, set the reaction's flux as a new objective, and ask the solver to maximize (resp. minimize) the flux for this reaction to get its upper (resp. lower) bound

This is what we call a flux variability analysis (see https://cobrapy.readthedocs.io/en/latest/simulating.html#Running-FVA). Note that, usually, at step 2 you would not constrain the model's previous objective to its optimal value, but to a fraction of its optimal in order to get the fluxes variability at, say, ±5% of the optimal growth rate.

Flux variability analysis can be done easily in cobra with (fraction_of_optimum=1.0 means that we keep 100% of the FBA's optimal objective value) :

solution_fva = cobra.flux_analysis.flux_variability_analysis(cmodel, fraction_of_optimum=1.0)

You can display unconstrained fluxes with something along solution_fva.loc[(solution_fva.maximum - solution_fva.minimum) > 0.01].

In your example, you will see that solution_fva.loc["LDH_Lm"] will return a minimal flux of -1000 and a maximal flux of 1000. Sooo it is clearly unconstrained: CPLEX "choose" a value of 0.16, but any value between -1000, 1000 is correct. Somehow, the .summary() change the value that is reported... but both are equally valid.

---

In terms of things we could do to improve COBRApy:

• maybe we could add (even more) explanations to the FVA page for newcomers to constraint-based models?
• maybe Reaction.flux should be private to avoid the temptation of using it?

HTH

[Midnighter]

As to why model.optimize() seemingly return the same flux distribution and do not stochastically select one value or another in the range of what is a solution...

Most solvers provide a hot start feature, i.e., a previous solution is used as the starting point for a new solution. Additionally, the same solver running on the same compute node, due to its implementation, will likely always arrive at the same solution. So yes, I highly recommend working with pFBA and FVA.

[Midnighter]

Other ways to get a sense of the non-uniqueness of flux distributions is to use geometric FBA or flux sampling.

[gcxxx]

@Midnighter Sorry to bother you again... I have thoroughly read the flux sampling section in the document, but I am not fully getting the sense of ACHR sample method. In one of the article I am reading, the author also use ACHR sample method,

But the sentence painted in yellow is confusing, how to sample 100000 points with 10000 iteration? thinning factor=0.1? What's the meaning of unique here?

[Midnighter]

Since our hope is that these questions on the discussion board will serve as a collection of knowledge for future users, can you please post this question in a new topic separately?

[gcxxx]

OK, since I am now focusing on another project, I will repost this question in the near future😂