Optimization error with VMH diets

[ryan4martin]

Hello,

I've been trying to run simulations with the diets other than the western diet (which runs without any problems) to work with MICOM and the AGORA models. To start I have been focusing on the VMH high fat low carbs diet you have provided that is compatible with MICOM.

The problem arises when I use the grow function with the VMH high fat low carbs diet. The majority of samples no longer are able to grow and I get the following error:

"Could not solve cooperative tradeoff for sample_x. This can often be fixed by enabling presolve, choosing more permissive atol and rtol arguments, or by checking that medium fluxes are > atol."

As the error suggestions, I tried the same setup with presolve enabled and the same error was returned.

I then found the following issue (micom-dev/paper#5) where you suggest directly running the cooperative_tradeoff function to get more specific error messages.

The following is the output I got from trying cooperative_tradeoff directly:

from micom import load_pickle
from micom.qiime_formats import load_qiime_medium

com = load_pickle("sample_x_model.pickle")
medium = load_qiime_medium("vmh_high_fat_low_carb_agora.qza")
com.medium = medium['flux']
sol = com.cooperative_tradeoff(fraction=0.5)

[19:38:52] WARNING  I could not find the following exchanges in your model:  community.py:668
                    EX_cellul_m, EX_datp_m, EX_nadp_m, EX_dtmp_m, EX_amet_m,                 
                    EX_tet_m, EX_starch1200_m, EX_cmp_m, EX_gtp_m,                           
                    EX_chsterol_m                                                            
[19:38:53] WARNING  solver encountered an error infeasible                    solution.py:220
           WARNING  solver encountered an error infeasible                    solution.py:220
OptimizationError: could not get community growth rate.
---------------------------------------------------------------------------
OptimizationError                         Traceback (most recent call last)
<ipython-input-27-c2178fd3e951> in <module>
      5 medium = load_qiime_medium("vmh_high_fat_low_carb_agora.qza")
      6 com.medium = medium['flux']
----> 7 sol = com.cooperative_tradeoff(fraction=0.5)

~/miniconda3/envs/qiime2/lib/python3.8/site-packages/micom/community.py in cooperative_tradeoff(self, min_growth, fraction, fluxes, pfba, atol, rtol)
    794             rtol = self.solver.configuration.tolerances.feasibility
    795 
--> 796         return cooperative_tradeoff(
    797             self, min_growth, fraction, fluxes, pfba, atol, rtol
    798         )

~/miniconda3/envs/qiime2/lib/python3.8/site-packages/micom/problems.py in cooperative_tradeoff(community, min_growth, fraction, fluxes, pfba, atol, rtol)
     77         com.objective = com.scale * com.variables.community_objective
     78         min_growth = (
---> 79             optimize_with_retry(com, message="could not get community growth rate.")
     80             / com.scale
     81         )

~/miniconda3/envs/qiime2/lib/python3.8/site-packages/micom/solution.py in optimize_with_retry(com, message)
    244         sol = com.optimize()
    245     if sol is None:
--> 246         raise OptimizationError(message)
    247     else:
    248         return sol.objective_value

OptimizationError: could not get community growth rate.

I am not sure what else to try in order to successfully implement the VMH diets. Any suggestions are greatly appreciated!

Thanks!!

[cdiener]

Hmm, I just tried with a model I had laying around and it did work:

In [1]: import micom as mm

In [2]: com = mm.load_pickle("A1.pickle")

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

In [4]: med
Out[4]: 
              reaction metabolite    global_id      flux
reaction                                                
EX_ala_L_m  EX_ala_L_m    ala_L_m  EX_ala_L(e)  3.519601
EX_arg_L_m  EX_arg_L_m    arg_L_m  EX_arg_L(e)  1.873749
EX_asp_L_m  EX_asp_L_m    asp_L_m  EX_asp_L(e)  4.159954
EX_btn_m      EX_btn_m      btn_m    EX_btn(e)  0.000100
EX_ca2_m      EX_ca2_m      ca2_m    EX_ca2(e)  1.677521
...                ...        ...          ...       ...
EX_gsn_m      EX_gsn_m      gsn_m    EX_gsn(e)  0.000148
EX_nadp_m    EX_nadp_m     nadp_m   EX_nadp(e)  0.000003
EX_urea_m    EX_urea_m     urea_m   EX_urea(e)  0.084613
EX_dgsn_m    EX_dgsn_m     dgsn_m   EX_dgsn(e)  0.000167
EX_nmn_m      EX_nmn_m      nmn_m    EX_nmn(e)  0.000529

[122 rows x 4 columns]

In [5]: com.medium = med.flux
[12:52:00] WARNING  I could not find the following exchanges in your model: EX_cmp_m, EX_amet_m                                                               community.py:668

In [6]: sol = com.cooperative_tradeoff(fraction=0.5)

In [7]: sol
Out[7]: <CommunitySolution 0.003 at 0x7fa6d5af9130>

The following info would be helpful:

• Are you using the genus-level for your models?
• Did you download the diet recently?
• What solver are you using?

[ryan4martin]

1. I am using strain level AGORA models.
2. The most recent time I downloaded the diet was today.
3. CPLEX

Update:
I tried with a genus-level model for the same sample and am able to successful simulate the community. I also tried the complete_db_medium to make sure all the strain-level agora models could grow, but ran into the same error with cooperative_tradeoff.

[cdiener]

Yeah, the media are only completed on the genus level. Probably better to do it on the strain level. I will switch to this strategy. complete_db_medium should work with the strain level. What output do you get from complete_db_medium with the strain models? Especially, manifest.can_grow.value_counts()? If you get taxa without growth there you might have to increase the basic fluxes in the template medium before filling in gaps.

This medium is definitely complicated though since it simply does not contain enough carbon sources for Bacteria to grow. I'm also planning to add in mucin to alleviate this a bit.

[ryan4martin]

manifest.can_grow.value_counts() True 815 False 3 Name: can_grow, dtype: int64

[cdiener]

That doesn't look too bad though. Can you paste the code you are using to complete the strain medium and use it?

[phylatechnologies]

I am starting with the current qiime medium available on this repo.

`
from micom.qiime_formats import load_qiime_medium
medium = load_qiime_medium('vmh_high_fat_low_carb_agora.qza')

from micom.workflows.db_media import complete_db_medium

manifest, imports = complete_db_medium("AGORA-with-MUCINS.zip", medium, growth=1e-4,
threads=59, max_added_import=10, weights="C")

import pandas as pd

fluxes = imports.max()
fluxes = fluxes[(fluxes > 1e-6) | fluxes.index.isin(medium.reaction)]
completed = pd.DataFrame({
"reaction": fluxes.index,
"metabolite": fluxes.index.str.replace("^EX_", "", regex=True),
"global_id": fluxes.index.str.replace("_m$", "(e)", regex=True),
"flux": fluxes
})

from qiime2 import Artifact

arti = Artifact.import_data("MicomMedium[Global]", completed)
arti.save("vmh_high_fat_low_carb_strain.qza")

To grow the samples:

medium = load_qiime_medium('vmh_high_fat_low_carb_strain.qza')
growth_results = grow(manifest=manifest, model_folder='sample_models', medium=medium, threads=59, tradeoff=0.5)
`

[cdiener]

I see. It looks like "AGORA-with-MUCINS.zip" is the AGORA file downloaded from VMH. This is by itself not a valid MICOM model database so I'm pretty sure the problem lies there. MICOM model databases are the ones that you can download from https://zenodo.org/record/3755182 or one that was created using build_database. It's definitely a bug in MICOM though since it should have detected that and complained.

The completed medium looks okay though. Hard to say what the problem is without seeing the actual models. This should only happen if the strains in your community are exactly the strains that could not grow in the optimization. How many strains are there for each model? Numerical issues usually get worse the larger the model gets, though we have run a few models with 100+ taxa and it usually worked...

[ryan4martin]

The "AGORA-with-MUCINS.zip" was built with MICOM's build_database.

The samples have on average 55 strains (max of 109) in each model. It doesn't look like the number of strains matched in the sample correlates with whether the optimization fails or not. Is there a way I could send a few samples that the optimization is failing? I am unable to attach them here.

[cdiener]

Ah, that makes more sense. Definitely, you can send me the file or a Google drive link to mail (a) cdiener.com. A pickle for a single sample that fails is probably enough.

[cdiener]

Thanks for the models. I checked them and surprisingly some of the AGORA models have demand reactions and since those can't be fulfilled in low nutrient conditions you get the infeasibility. So no numerical issues, the models were indeed infeasible 😅 For instance, for sample 1, the following fixes the issue:

com.reactions.DM_btn__Faecalibacterium_prausnitzii_M21_2.lower_bound = 0
com.reactions.get_by_id("DM_thmpp(c)__Lactobacillus_fermentum_ATCC_14931").lower_bound = 0

This should be addressed in MICOM by adjusting those bounds. Having growth rate-independent demands is not good practice, This can also trip up complete_medium, so this should also work better after that. Hopefully, I'll find some time to get to this soon.

[ryan4martin]

Thank you! I am now able to simulate all my samples by looking for reactions with positive lower bounds or negative upper bounds and changing the bounds to 0.

[cdiener]

Awesome, looks like this was the issue then.

Now tracked in micom-dev/micom#56