Different outputs using different solver

[vgenisel]

Hello,
I've been playing with MICOM and noticed that simulations using the same conditions (taxa abundance and diet) but with different solvers generate different growth rates and fluxes of produced metabolites. I understand that small changes are expected considering the different accuracies of the solvers but in some cases even the profile of the fastest-growing bacteria and the most-produced metabolites changes.
Does anyone observe the same?
I'm using CPLEX 22.1, GUROBI 10.0.1 and the default OSQP solver.
Cheers

[cdiener]

Hi, in general I would say yes, we have observed something similar, but some part of that is expected. Can you give some detail how you got the fluxes in the end? If you used the grow workflow what strategy did you specify? In general flux distributions are not unique given a single unique optimum, so different flux distributions can give you the same optimal solution. Strategies like minimal medium and pFBA help to reduce that variability but make no guarantees.

Having said that, the growth rates have a unique solution. So if those vary that is probably a numerical issue. OSQP struggles a lot with numerically difficult problems, especially linear ones. Playing around with the iteration limit can help there but I have not found a great default there. Some problem improve with more iterations but others get suber bad. CPLEX and Gurobi are generally more accurate there. If the growth rates between CPLEX and Gurobi are very different that is probably something we want to try to address/fix. Looks like your data is publicly available, so would you mind sharing the taxonomy/abundance tables (the input to grow or Community)?

[vgenisel]

Thanks for the reply.
I used the grow workflow with the default condition (pFBA strategy). The tradeoff was set to 0.5 for the 3 solvers. Then I used the function production_rates to generate the fluxes.

I noticed that GUROBI and CPLEX agree better for the fluxes of metabolites but are not consistent with each other for the growth rates. Some taxa grow much more using GUROBI than CPLEX (around 10 times), while others have the growth compromised due to the dominant bacteria and don't even grow at all (the same happens in CPLEX too). This made me curious, considering that the medium was designed using the function complete_db_medium to ensure a growth rate = 0.01. I believe it's a numerical issue, indeed.

Concerning the data, I used a VMH high_fibre diet, taxonomy from the ISB 2020 workshop (only for the 4 healthy individuals but I can't attach qiime artifacts) and AGORA reconstructions on genus level (qiime artifact, found then in the MICOM repository https://zenodo.org/record/3755182#.Y_PgGSZByUk).

I ran other simulations using taxonomy from infants at weaning age and the same behaviour with the solvers was observed (results in the ppt attached).

high_fibre.csv
Gurobi vs OSQP vs CPLEX MICOM infants 100223.pptx

[cdiener]

Perfect, thanks, I will investigate. What version of MICOM did you use?

[vgenisel]

Great, thanks. It was the version 0.32.3

[cdiener]

Hi, so I checked one of the models and I could reproduce your results. Checking the solver status Gurobi marked the problem as numerically unstable so it is likely due to numerical issues. It seems to stem from the fact that the high fiber diet you used allows for only very small growth rates. If I run it with the flux bounds from the medium multiplied by 100 all of the solvers converge with an "optimal" status. Alternatively you can also use one of the more recent media from the MICOM media repository, for instance when using the high fiber medium from there you get:

Python 3.11.2 (main, Feb 12 2023, 00:48:52) [GCC 12.2.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.10.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import micom as mm

In [2]: med = mm.qiime_formats.load_qiime_medium("/home/cdiener/Downloads/vmh_high_fiber_agora.qza")

In [3]: sols = {}

In [4]: for so in ["cplex", "gurobi", "osqp"]:
   ...:     com = mm.load_pickle(f"{so}.pickle")
   ...:     com.medium = med.flux
   ...:     sol = com.cooperative_tradeoff(fraction=0.5)
   ...:     print(com.taxonomy.sample_id[0], com.solver.status, flush=True)
   ...:     sols[so] = sol.members
   ...: 
[13:26:52] WARNING  I could not find the following exchanges in your model: EX_urate_m, EX_ttdcea_m,   community.py:669
                    EX_hdcea_m, EX_10fthf_m, EX_ascb_L_m, EX_caro_m, EX_arach_m, EX_xylt_m, EX_octa_m,                 
                    EX_retinol_m, EX_phyQ_m                                                                            
ERR1883210 optimal
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 937956
Academic license - for non-commercial use only - registered to [email protected]
Read LP format model from file /tmp/tmpns57yi7a.lp
Reading time = 0.28 seconds
: 41704 rows, 124262 columns, 516998 nonzeros
[13:27:15] WARNING  I could not find the following exchanges in your model: EX_urate_m, EX_ttdcea_m,   community.py:669
                    EX_hdcea_m, EX_10fthf_m, EX_ascb_L_m, EX_caro_m, EX_arach_m, EX_xylt_m, EX_octa_m,                 
                    EX_retinol_m, EX_phyQ_m                                                                            
ERR1883210 optimal
[13:27:37] WARNING  I could not find the following exchanges in your model: EX_urate_m, EX_ttdcea_m,   community.py:669
                    EX_hdcea_m, EX_10fthf_m, EX_ascb_L_m, EX_caro_m, EX_arach_m, EX_xylt_m, EX_octa_m,                 
                    EX_retinol_m, EX_phyQ_m                                                                            
ERR1883210 optimal

In [5]: for so in sols:
   ...:     obj = sols[so].growth_rate.pow(2).sum()
   ...:     diff = (sols[so].growth_rate - sols["cplex"].growth_rate).abs().mean()
   ...:     print(f"solver: {so} - objective: {obj} - mean diff to ref solution: {diff}")
   ...: 
solver: cplex - objective: 0.03034629534902814 - mean diff to ref solution: 0.0
solver: gurobi - objective: 0.029833997946152057 - mean diff to ref solution: 0.003064592516437963
solver: osqp - objective: 1.1368086781097095 - mean diff to ref solution: 0.04743052249275915

OSQP still does not do a good job but the other solvers mostly agree with each other. If your medium is very limited you could also multiply the initial fluxes with some larger number (10-100) but would then have to adjust your output to get an approximation for the original medium (dividing by that same scale). The media in the repository are usually vetted to allow for sufficient growth so they should mostly work.

[vgenisel]

Hi,
Many thanks for finding a solution. It is good to know that the scale of the dietary fluxes may influence the predictions.
When designing my medium, I imported the fluxes from the VMH database but divided them by 24 (to convert mmol/d to mmol/h). Probably that was the origin of the issue.