minimal_medium() returns more reactions than initially present

Dear maintainers,

I'm encountering a strange issue with micom.media.minimal_medium() in MICOM (version 0.16.1). Apparently the medium computed contains more metabolites than the initial medium provided to the community. From what I understand, this should not be possible with the parameter open_exchanges=False (more details below).

I have a community of roughly ~250 species from AGORA 1.03 (without mucins) and a custom medium derived from the WesternDiet provided in the original MICOM publication: western_diet_custom.txt. I also use CPLEX 12.9.0.0.

I apply the medium on my community in a similar way than shown in the documentation:

# List of exchange reaction in the model
ex_ids = [r.id for r in com.exchanges]
# Find medium reaction that are recognized by the community
medium_in_ex_ids = medium.index.isin(ex_ids)
# Define medium of community
com.medium = medium[medium_in_ex_ids]

Then, I run a cooperative tradeoff to obtain the community growth rate:

# Compute growth rate on diet
sol = com.cooperative_tradeoff(fraction=0.7, fluxes=True, pfba=False,)
sol

So far so good. And I'm also able to compute the minimal medium necessary to obtain 95% of this community growth rate:

# Compute minimal medium necessary to reach 95% growth rate
med = micom.media.minimal_medium(
    com,
    0.95 * sol.growth_rate,
    exports=False,
    open_exchanges=False,  # !!!
)
med

EX_26dap_M_m     0.013413
EX_2dmmq8_m      0.011625
EX_acgam_m      23.994033
EX_adocbl_m      0.000054
EX_arg_L_m       3.600214
                  ...    
EX_isochol_m     0.000001
EX_HC02194_m     0.000033
EX_uchol_m       0.000032
EX_sel_m         0.001193
EX_cellul_m      0.000836
Length: 130, dtype: float64

This is where I notice a problem: when I compare the fluxes of the minimal medium with the original medium, there are a handful of extra reactions that were not initially present (reproduced below). From what I understand, it should not be possible with open_exchanges=False. Any idea of what is happening? I hope I'm giving enough information (feel free to ask for more, I didn't share the pickle of the community because it's a 300MB file, but I could find a way if it's necessary).

Note that I did not observe this problem with AGORA 1.03 (with mucins) or AGORA 1.02, and it seems to only happen with some specific samples. My main goal was to somehow use minimal_medium() in order to identify essential metabolites required by my communities, but the fact that the function add more reactions than expected might be troublesome.

Cheers, Nils

This has to do with tolerances in the solver. A numerical solver is bound by floating point accuracy, thus it will never return the exact "true" solution of the problem but rather a solution within a certain tolerance. This is usually managed by a relative and absolute tolerance. What you are seeing here is the effect of the relative tolerance eps_rel. The solver usually will only guarantee that your solutions will be within eps_rel * ||v||* where ||.||* is the infinity norm (maximum absolute value). The default value for eps_rel with cplex is 1e-6 and you have at least one flux with a value of around 20, so your solution will be within a 1e-4 to 1e-6 of the real solution. All of those impossible imports are thus not distinguishable from zero for the solver.

I am working on an updated version that lets you define the threshold for returning imports. For now minimal_media returns any import flux > 1e-6 but you will be able to change this cutoff soon. To alleviate this you can either drop imports below a certain threshold or decrease the solver tolerance with com.tolerance = 1e-8, for instance. However, this will affect performance and may affect convergence. To put it into perspective, single bacteria have fL volumes and the default flux unit of micom is mmol/h. So that corresponds to less than a single molecule imported per cell and hour. This is unlikely to have a large impact on downstream analyses.

P.S.: I would also set min_growth=0.95*sol.members.growth_rate.drop("medium") in minimal_media. Otherwise you are loosing all the advantages from the cooperative tradeoff.

Hi @cdiener,

Thank you for your answer, that makes perfect sense. Indeed I have very large values for my fluxes (hence very large error tolerance), this is because I don't use mmol/gDW/h but rather mmol/human/day. The reasoning was that, since AGORA models do not have any internal constraints, this should be transparent for the community model (the outputs should just be in the same units than the inputs, with no unintended consequences on the flux distribution). Since the goal was to compare results with other community modeling approaches that use mmol/human/day, but also to compare different literature diets that are reported in mmol/human/day (the VMH diets for instance), I wanted to avoid unnecessary conversions.

From your answer, I now realize though that MICOM might have some default thresholds that do expect fluxes in the range of mmol/gDW/h. Maybe I should change my workflow and actually convert my diets before running MICOM, and convert back my outputs after if I want to compare with other methods.

P.S.: I would also set min_growth=0.95*sol.members.growth_rate.drop("medium") in minimal_media. Otherwise you are loosing all the advantages from the cooperative tradeoff.

Indeed, I checked back the documentation and it's there. I think months ago I initially adapted what you did in the original MICOM paper https://github.com/micom-dev/paper/blob/master/workflows/media_and_gcs.py#L19-L41, which does not use the min_growth parameter. I guess this functionality has been added later and I didn't notice. I'm really grateful for the tip, and will adapt my scripts.

Thank you again for your help.

Nils

micom-dev / micom

minimal_medium() returns more reactions than initially present #15