Solution fluxes that were not in com.medium

[jmol0917]

Hi All,

I am concerned I might have a but in my code.

I am working with a 3-strain community in modified M9 media:

m9_media = [
    'EX_ca2(e)', 'EX_cl(e)', 'EX_cu2(e)', 'EX_na1(e)', 'EX_nh4(e)',
    'EX_k(e)', 'EX_pi(e)', 'EX_h2o(e)',
    'EX_mg2(e)', 'EX_so4(e)', 'EX_zn2(e)', 'EX_o2(e)', 
    'EX_mn2(e)', 'EX_cobalt2(e)',                            # These are not in the paper 
    'EX_fe3(e)', 'EX_fe2(e)'                                 # In small amounts in casamino acids
]                                                                        

casamino_acids = [
    'EX_glu_L(e)', 'EX_pro_L(e)', 
    'EX_ala_L(e)', 'EX_ser_L(e)', 'EX_phe_L(e)', 'EX_gln(e)', 'EX_thr_L(e)', 'EX_val_L(e)',
    'EX_met_D(e)', 'EX_met_L(e)', 'EX_leu_L(e)', 'EX_lys_L(e)', 'EX_ile_L(e)', 'EX_tyr_L(e)',
    'EX_gly(e)', 'EX_Lcyst(e)', 'EX_arg_L(e)', 'EX_asn_L(e)'
#     'EX_asp_L(e)', 'EX_his_L(e)'  These two are also missing according to some sources
    # No tryptophan in casamino acids
]


carbon_source = 'EX_glc_D(e)'

customised_m9 = {}

for nutrient in m9_media:
    customised_m9.update({nutrient.replace('(e)', '_m'): 1000})
for nutrient in casamino_acids:
    customised_m9.update({nutrient.replace('(e)', '_m'): 40})
customised_m9.update({carbon_source.replace('(e)', '_m'): 10})

Then I pass it to community.medium:

with community:
    medium = community.medium

    for reaction in medium:
        if reaction not in customised_m9:
            medium[reaction] = 0.0
        else:
            medium[reaction] = customised_m9[reaction]
    community.medium = medium
    community.add_cons_vars(constraint)
    solution = community.cooperative_tradeoff(fluxes=True, pfba=True)

Where there is no EX_pheme(e) for example, but when I check export fluxes all of the strains seem to be consuming it:

Is there something I am missing? Do negative signs actually indicate export? Because I checked and pheme biosynthesis within the three is much lower than the values indicated above.

Thanks!

[cdiener]

Hmm, I don't see anything wrong right now. Maybe double check your community.medium before the cooperative_tradeoff call? If you are using OSQP that can sometimes give very inaccurate solutions, though this would be a pretty large error. What does ,solution.status look like? Also make sure you don't run another cooperative_tradeoff outside of the with: block as this will reset back to the previous medium.

[jmol0917]

Thanks for your reply @cdiener.

I am using OSPQ: <optlang.osqp_interface.Model at 0x1fbf98c9f70>, I guess that is the problem? solution.status is iteration_limit. I am not running other cooperative_tradeoff outside the with: block.

I am checking community.medium by printing it right before running community.cooperative_tradeoff. Is there a better way of double checking?

Thanks!

[cdiener]

Printing it right before should be enough. Or just checking community.reactions.EX_pheme_m.bounds. If those are indeed zero it is the solver. You could try it with setting community.solver.configuration.presolve = True and see if it gets better.

If you are at an educational institution you could also try with a free academic license for cplex or gurobi...

[jmol0917]

Hmmm, when I do community.reactions.EX_pheme_m.bounds I get (-0.0, 1000.0) when checking within with com: and (-1000.0, 1000.0). I have checked with a few other reactions and they are all the same for those I don't want in my media.

However, when I manually change EX_pheme_m bounds to be (-0.0, 0.0) I get very similar results in terms of fluxes. So I guess it is the solver in the end?

Also, do I need both bounds to be (-0.0, 0.0)? It seems like the with com: manipulation of com.medium only changes the lower bound to produce (-0.0, 1000.0)

I am in talks with Gurobi to get one of their licenses.

Thanks.

[cdiener]

I should mention that finding a better open source solver is definitely on my list. For small models like yours it would also be feasible to do a crossover with a different LP solver like GLPK to get a solution with better feasibility. The trick would be to take the growth rate estimates from the cooperative tradeoff and then limit the growth rates to those and do the pFBA with a different solver.

[jmol0917]

I got a Gurobi licence and that addressed the issue I was having here. Thanks!