medium in grow/growth vs medium in single model workflow

[jmol0917]

Hi!

I ran a few samples using your multi-sample workflow (build/grow) and a defined medium and everything runs with no issues - I get my predictions and all. Then, I tried grabbing the models from one of my samples and run things based on Single Model - build the exact same community, with the same abundances and same medium. However, during the later I keep getting an OptimizationError: could not get community growth rate

This only happens when I pass the same medium I used for my multi-sample workflow. If I leave all the exchanges open, I do get growth predictions with Community.cooperative_tradeoff. In both cases I am using the same tradeoff value of 0.7

I guess this is happening because there are essential nutrients missing from my medium. My question is why are things working just fine when I run the multi-sample workflow? Is media processed differently or supplemented at one point in the workflow? Apologies if I am missing something. I did have a look at the source code, and as far as I understand the media is not modified.

I am using the latest version of Gurobi as my solver.

Many thanks!
JP

[jmol0917]

For reference, this is how I am passing my medium in each case:
For multi-sample modelling:

medium_df = pd.read_csv('medium_file.csv', index_col = 0)
samples = micom.workflows.grow(manifest,
                              'path/to/models/',
                              medium_df, 
                              tradeoff = 0.7,
                              strategy = 'none',
                              threads = 10)

Single sample modelling:

medium_df = pd.read_csv('medium_file.csv', index_col = 0)
medium = dict(zip(medium_df['reaction'], medium_df['flux']))
comm.medium = medium
comm.solver.configuration.presolve = True
coo_solution = comm.cooperative_tradeoff(0.7, fluxes=True, pfba=True)

Thanks!

[cdiener]

Seems like something is different in the building. If you use the model from the build workflow does it work? Like:

from micom import load_pickle

comm = load_pickle("path/to/models/A1.pickle")
comm.medium = med.flux
[...]

[jmol0917]

Just gave that a try and I am getting the same error. I tried different samples and the same thing is happening.

This is the code I am currently using. I am processing the media differently because process_medium was asking for a second argument (sample id) and my media is not formatted that way:

import micom
import pandas as pd

comm = micom.load_pickle("models/studyA/D1.pickle")      # Using a different community this time
medium = pd.read_csv('studyAmedium.csv', index_col = 0)
inulinMedium = pd.concat([medium,
                           pd.DataFrame({'reaction':['EX_inulin_m'],
                                         'metabolite':['inulin_m'],
                                         'global_id':['EX_inulin(e)'],
                                         'flux':[14]
                                        },index = ['EX_inulin_m'])])

inulin_modelling_medium = dict(zip(inulinMedium['reaction'], inulinMedium['flux']))

with comm:
    print(len(comm.medium))
    medium = comm.medium
    for element in medium:
        if element in inulin_modelling_medium:
            medium[element] = inulin_modelling_medium[element]
         else:
            medium[element] = 0.0
  
    comm.medium = medium
    print(len(comm.medium))
    comm.solver.configuration.presolve = False
    coo_solution = comm.cooperative_tradeoff(0.7, fluxes=True, pfba=True)

If I run coo_solution = comm.cooperative_tradeoff(0.7, fluxes=True, pfba=True) without specifying the media things run just fine.
I was going to upload the pickle file I am using but the format is not compatible with GitHub.

Thanks for any help!

[cdiener]

Hi and so sorry for the delay. I did run your example files for sample A1 and everything seems to work for me. Here is what I did:

In [13]: import micom as mm

In [14]: import pandas as pd

In [15]: ab = pd.read_csv("example_complete_taxa_updatedPath.csv", index_col=0).drop(columns="file")

In [16]: ab
Out[16]: 
                 sample_id           species           phylum  ...               genus abundance                id
equolifaciens           A1     equolifaciens   Actinobacteria  ...       Adlercreutzia     0.011     equolifaciens
muciniphila             A1       muciniphila  Verrucomicrobia  ...         Akkermansia     0.024       muciniphila
hadrus                  A1            hadrus       Firmicutes  ...        Anaerostipes     0.039            hadrus
[...]

[32 rows x 9 columns]

In [17]: com = mm.Community(ab, model_db="model_db")
[09/17/24 09:29:48] WARNING  The following sinks or demands were forced in the model and have been relaxed since they would    util.py:64
                             interfere with MICOM's growth rate estimation: DM_btn                                                       
Building ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:01:04

In [18]: medium = pd.read_csv('studyAmedium.csv', index_col = 0)
    ...: inulinMedium = pd.concat([medium,
    ...:                           pd.DataFrame({'reaction':['EX_inulin_m'],
    ...:                                          'metabolite':['inulin_m'],
    ...:                                          'global_id':['EX_inulin(e)'],
    ...:                                          'flux':[14]
    ...:                                        },index = ['EX_inulin_m'])])

In [19]: com.medium = inulinMedium.flux

In [20]: sol = com.cooperative_tradeoff(fraction=0.7, fluxes=True, pfba=True)

In [21]: sol.members
Out[21]: 
                  abundance  growth_rate  reactions  metabolites
compartments                                                    
adolescentis       0.029058     0.045340       1129          991
aerofaciens        0.060120     0.092943        920          850
bicirculans        0.087174     0.006905       1032          926
bolteae            0.006012     0.016527       1932         1450
cellulosilyticus   0.011022     0.017973       2278         1513
champanellensis    0.003006     0.058372       1404         1140
cylindroides       0.008016     0.000000        994          892
difficile          0.028056     0.043676       1844         1401
dorei              0.023046     0.035818       2345         1467
equolifaciens      0.011022     0.022707        942          904
faecium            0.006012     0.040897       1528         1175
fragilis           0.036072     0.056031       2500         1533
gallolyticus       0.002004     0.000000       1225         1076
gnavus             0.012024     0.021742       2225         1454
hadrus             0.039078     0.000000       1704         1331
hansenii           0.006012     0.028931       2055         1297
longum             0.004008     0.012719       2648         1499
medium                  NaN          NaN        758          758
muciniphila        0.024048     0.037401       2274         1386
ovatus             0.028056     0.043509       2673         1550
piger              0.010020     0.023326        980          890
plautii            0.005010     0.000000       1207         1105
prausnitzii        0.163327     0.252244       1972         1467
producta           0.003006     0.019832       1235         1092
scindens           0.003006     0.026169       2049         1394
smithii            0.028056     0.043513        714          742
sporogenes         0.014028     0.022307       1766         1356
stadtmanae         0.003006     3.826743        823          809
thetaiotaomicron   0.020040     0.031754       2524         1549
uniformis          0.126253     0.195106       2418         1500
vulgatus           0.165331     0.256239       2474         1519
xylaniphila        0.010020     0.019706       1930         1356
xylanisolvens      0.025050     0.038884       2322         1422

In [23]: sol.fluxes.loc[:, "EX_but(e)"]
Out[23]: 
compartment
adolescentis         -0.007000
aerofaciens                NaN
bicirculans                NaN
bolteae              -0.000073
cellulosilyticus           NaN
champanellensis            NaN
cylindroides        256.351269
difficile           356.287920
dorei                      NaN
equolifaciens              NaN
faecium                    NaN
fragilis                   NaN
gallolyticus               NaN
gnavus              -56.454409
hadrus              204.891429
hansenii                   NaN
longum                     NaN
medium                     NaN
muciniphila                NaN
ovatus                     NaN
piger                      NaN
plautii             350.968587
prausnitzii         241.337914
producta                   NaN
scindens                   NaN
smithii                    NaN
sporogenes          271.039608
stadtmanae                 NaN
thetaiotaomicron           NaN
uniformis                  NaN
vulgatus                   NaN
xylaniphila                NaN
xylanisolvens              NaN
Name: EX_but(e), dtype: float64

The file column isn't needed as long as you use a model db even if you unzipped it. My model_db folder is the data folder from the agora103_species.qza after unzipping. Similar for the medium setting it is enough to set with any dict-like object, like a pandas Series or DataFrame column. MICOM will automatically figure out if reactions from your medium are missing and set the ones that were found.

Not sure where the discrepancy is coming from.

[jmol0917]

Thanks heaps for looking into it!

I tried replicating your code step by step and I still got an error - this time a GurobiError (which, in my head, is not too different from the error I was getting before), so I think this has more to do with the solver, which is still interesting because things do work when using the build/grow workflow.

Just in case, the version of gurobi I am working with is 11.0.0

[cdiener]

Is the error you are seeing something along the line of cant access "X" on Model. Then I can reproduce it. Seems to be numerical issues which occur only when using pfba=True as this won't attempt a crossover from the QP to LP first like the other approaches. It should work right away with pfba=False. It also works for me when setting com.solver.problem.Params.NumericFocus = 2 before the cooperative tradeoff call. I will send in a fix though since that also affects the hybrid solver where it even returns a wrong solution.