Varying biomass trends with different sequence of model additions

[SandhyaVasu]

Hello @dukovski @liamkelley93 @jeremymchacon @mquintin @dileep-kishore

Firstly, I would like to begin by saying - COMETS is just awesome! I am also very happy with your team for the earlier prompt response.
I am aware that a few years back, a similar issue was raised on this platform, but I wanted to clarify with you guys again (as you had suggested to go with a small timestep, but it is not helping).

I am performing microbial community simulations wherein I want to mimic a condition where 3 (or more) microbes enter an environment simultaneously.
However, the order of the models in the list : [1, 2, 3] ; [2, 1, 3] produces drastic differences in the growth trends for certain permutations (please check below). These are for a very small timestep - 0.01. (Attached at the end are the parameters I used.) I am using the Agora high fibre diet. It is in this format: {'arab_L': 0.04736842}

Note: I am adding all models at t=0 as follows:

model_list = [1,2,3] #assume 1,2,3 are comets models
comm_layout = c.layout()
for i in model_list:
comm_layout.add_model(i)

My current goal is to add models to a community layout and test if all members survive or not. But since this is stochastic (as in, order decides the outcomes), I am trying to perform all possible sequences of addition for a given set of three members and trying to compute the probability of all species surviving. But this leads to performing n! simulations for every set, and COMETS crashes when I run back-to-back simulations in a loop.

Can you please assist me in using COMETS in such a way that the order does not matter? This will be really, really helpful! COMETS is the only best dFBA tool I have come across for modelling microbial communities and if you can help me solve this issue, it will benefit me a lot!

Thanks and Regards
Sandhya

These are the parameters I set:

comm_params.set_param('randomSeed', 42)
comm_params.set_param("spaceWidth", 20)
comm_params.set_param("defaultVmax", 20.0)
comm_params.set_param("defaultKm", 5e-5)
comm_params.set_param('timeStep', 0.01)
#comm_params.set_param('maxSpaceBiomass', 10.0)
comm_params.set_param('maxCycles', 10)
#comm_params.set_param('growthDiffRate', 1)
comm_params.set_param('writeFluxLog', True)
comm_params.set_param('writeMediaLog', True)
comm_params.set_param('FluxLogRate', 1)
comm_params.set_param('MediaLogRate', 1)
comm_params.set_param('writeTotalBiomassLog', True)
#comm_params.set_param("minSpaceBiomass", 1e-3) #gDW
#Set chemostat dilution rate
dilution_rate = 0.017
comm_params.set_param('deathRate', dilution_rate)
comm_params.set_param('metaboliteDilutionRate', dilution_rate)

[dukovski]

Hi Sandhya,

Thank you for the kind words. This behavior is not supposed to happen. I think we fixed this bug, but it seems to be showing up again. Can you please send me the script and models so I could reproduce the bug and work on it?

Thanks,
Ilija

[SandhyaVasu]

Hi Ilija,

Thanks for getting back so soon!

The Jupyter notebook (containing the script) , the diet files, model files, and the YAML environment file are all attached below here as well as in this Google Drive folder: Link to files

The cometspy version I am using is: 0.5.2

Jupyter notebook with the code:
comets_model_additions_permutations.ipynb

Model files:

Lactobacillus_brevis_DSM_20054.xml
Faecalibacterium_prausnitzii_ERR1022279.xml
Bifidobacterium_adolescentis_ATCC_15703.xml

Diet file:

Magnusdottir_2017_TableS12_diets.csv

Conda environment - YAML file (exported from a Mac OS)
comets.yml

Thanks and Best Regards
Sandhya

[dukovski]

Hi Sandhya,

I am working on your problem. I am trying batch culture, I turned off the chemostat setting for simplicity. Here's a notebook with all metabolites set to 100. The runs are exactly the same.

comets_model_additions_permutations.ipynb

[dukovski]

If I set the initial values for all metabolites to a small value, and then set the others to their actual value, same thing, all good.

comets_model_additions_permutations.ipynb

[dukovski]

So you have a system here where some of the organisms cannot live without a nutrient produced by other. Now, if you order the organism that produces the limiting nutrient first, the second one will grow. But if you order so that the first one is the one that cannot live without the produced nutrient, it will not grow, and not produce the nutrient for the third etc. The solution is to kickstart your system with a tiny amount of all nutrients. This is actual y how should be done experimentally too. Here's my notebook to do that. Just add this code before adding the nutrients and refresh rates:
for met in comm_layout.media['metabolite']:
comm_layout.set_specific_metabolite(met, 0.001)

Here's the notebook, check cell 10.

comets_model_additions_permutations.ipynb

[SandhyaVasu]

Awesome, it now works!!! 👍 👍
Thank you so much, Ilija, for kindly assisting with this issue and taking the time to troubleshoot. It really means a lot!!

Please check if my understanding is right.

1. Despite adding the three models together at t=0, there is a small bias with respect to who arrives first. If the one who comes first is dependent on the second's metabolite(s), it fails to grow. So, to compensate for this, we add a tiny amount of all nutrients in the beginning.
2. Despite this, again, the first-arriver may utilise the initial resources in its own way, which results in different biomasses of each microbe across the different orders of addition. But all microbes grow, and the extent to which they grow depends on who came first and how they used the initial resources.

3. But, if we add a small amount of all metabolites at the start, will it not always result in all the microbes growing, as we are mimicking something like a minimal medium custom-made for all the models? So, will I be able to answer questions like: "Given a specific diet, how many microbes from A,B,C can successfully grow when I inoculate them at the same time in a gut-like environment?"
   In other words, will I be able to catch communities where some microbes don't grow at all?

4. I tried what you suggested with two models out of the three I had shared.

i) Setting the initial values for all metabolites to a small value: Consistent trend -> both grow, but to different extents

ii) Without setting the initial values for all metabolites to a small value : Contrasting trends

Thanks again!
Sandhya

[dukovski]

Yes to 1). In 2) the results should be the same regardless of the order. 3) Yes you will be able to catch that if you initialize with a small amount. They may grow at first a tiny bit, but will stop when the nutrient is gone. Do not replenish in a chemostat with everything, only with the ones in the media. Only start with tiny amount of all exchanged to kickstart the process.
4) In this case a) you either have too little initial media, or your time step is too large. The idea is that one strain will take some media to start growing, but if the time step is too large it one will take all of the initial, or not, depending on the order. Try putting large amount for the media, then go down with it until you start to see a difference.

[SandhyaVasu]

Okay, thank you for the clarification!
Will experiment with these and get back if I need help.

Regards
Sandhya

[SandhyaVasu]

Hi Ilija,

Hope you are doing well. Firstly, would like to thank you again for COMETS while beginning this post. Was re-visiting this problem that I had encountered before and it seems to be stemming from the below sorting issue in the Java files. Please correct me if I am wrong.

Quick summary of the issue

When two (or more) species are added to the same layout at t=0, simply changing the order
of the models in the input array (e.g. [A, B] vs [B, A]) produces different simulation
trajectories, even though the community composition is identical. We traced this to an
internal inconsistency between how medium metabolite names and medium refresh rates
are indexed.

Root cause

In FBAWorld.java, when the world's metabolite list is (re)built from the constituent
models, the metabolite names are explicitly alphabetically sorted:

// FBAWorld.java
Arrays.sort(newMetabNames);

However, the refresh rates (and related per-metabolite arrays supplied via the layout,
e.g. media_refresh) are populated in the order the metabolites were originally
encountered/added — i.e. the order in which models were passed in — rather than being
re-indexed to match the alphabetically sorted metabolite names.

The practical effect: refresh rate i is assumed to correspond to mediaNames[i], but
after the alphabetical sort, mediaNames[i] may no longer be the metabolite the refresh
rate was originally specified for. Since the input model array order changes the order in
which new metabolites are first encountered, it also changes which refresh rate ends up
misaligned with which metabolite — producing order-dependent results for what should be an
order-independent operation (adding the same set of models to a layout).

This affects FBAWorld.java and the layout-construction path in FBACometsLoader.java.

Test simulations

• Build a layout with two models A and B, non-trivial media_refresh rates.
• Run once with model array [A, B], once with [B, A].
• Compare biomass/media trajectories — they diverge, sometimes substantially (see attached
  before/after figure from a 20-random-pair sweep; several pairs show large qualitative
  differences in final community composition depending on array order alone).

Workaround we are currently using

On our Python script, sorting the media table alphabetically by metabolite name
before it is passed into comets() resolves the mismatch, since it forces the refresh
rates to already be in the same order the Java layer will sort them into:

layout.media = layout.media.sort_values('metabolite').reset_index(drop=True)

We also found it necessary to set maxSpaceBiomass to a high value (e.g. 1e40) to avoid
a separate biomass-cap interaction that was previously masking/compounding this issue in
our tests. With both changes, array order gives near-identical results, with residual
differences only at floating-point precision.

The attached figures show 20 randomly chosen species pairs simulated in both array orders ([A,B] solid, [B,A] dashed) over a 20h chemostat run. Before the fix, many pairs show visibly divergent biomass trajectories between the two orderings; after applying the sort fix, the solid and dashed curves for each pair overlap almost exactly.

Before sorting

After sorting

Happy to share our test script and the full before/after sweep data if helpful.

Thank you for maintaining COMETS!

Best,
Sandhya