MICOM Batch API

The basic unit of analysis in MICOM is called a “Batch”. It is ambiguous on purpose an can mean the following things:

• a single microbial community (like a batch of microbes, or a batch culture)

• a data set containing multiple microbial communities/samples (batch/group of samples)

Those are managed by the Batch class which contains all necessary tools to simulate those communities. It is the recommended entry point for MICOM as tries to strike a good balance between sensible defaults and flexibility in your analysis.

A batch follows this lifecycle:

1. Create a Batch from a taxonomy table, a growth medium, and either a model database or model files.

2. Call build() to assemble one community model per sample and write a manifest.

3. Call grow() to simulate growth, then inspect the returned GrowthResults.

4. Use tradeoff(), minimal_medium(), or complete_medium() for follow-up analyses.

Each batch comes with a coupled Configuration object which controls all details of the functions in the batch and that can be exported and saved to reproduce analyses in the future.

Like mentioned above a Batch requires the following:

1. A table of taxon/microbe abundances.

2. [optional] a growth medium/environment

3. [optional] a model database to link taxa to metabolic models

So let’s go through those.

Building community models

Taxon abundances

The first thing we need is a table of taxon abundances. As mentioned before this a table containing taxa and abundances. Best to look at an example first.

from micom.data import test_data

data = test_data()
data

	id	genus	species	strain	reactions	metabolites	abundance	sample_id
0	strain_a	Escherichia	Escherichia coli	strain a	95	72	764	sample_1
1	strain_b	Escherichia	Escherichia coli	strain b	95	72	783	sample_1
2	strain_c	Escherichia	Escherichia coli	strain c	95	72	154	sample_1
3	strain_d	Escherichia	Escherichia coli	strain d	95	72	98	sample_1
4	strain_a	Escherichia	Escherichia coli	strain a	95	72	86	sample_2
5	strain_b	Escherichia	Escherichia coli	strain b	95	72	212	sample_2
6	strain_c	Escherichia	Escherichia coli	strain c	95	72	518	sample_2
7	strain_d	Escherichia	Escherichia coli	strain d	95	72	261	sample_2
8	strain_a	Escherichia	Escherichia coli	strain a	95	72	12	sample_3
9	strain_b	Escherichia	Escherichia coli	strain b	95	72	170	sample_3
10	strain_c	Escherichia	Escherichia coli	strain c	95	72	51	sample_3
11	strain_d	Escherichia	Escherichia coli	strain d	95	72	799	sample_3
12	strain_a	Escherichia	Escherichia coli	strain a	95	72	878	sample_4
13	strain_b	Escherichia	Escherichia coli	strain b	95	72	996	sample_4
14	strain_c	Escherichia	Escherichia coli	strain c	95	72	156	sample_4
15	strain_d	Escherichia	Escherichia coli	strain d	95	72	150	sample_4

The table may contain additional taxonomy ranks. MICOM uses every supplied rank to make database matching more stringent. For example, supplying both genus and species requires a model to match both ranks. Omit a rank when that annotation is unreliable.

The required columns are id, sample_id, abundance, and the taxonomy rank used by the model database. In this example, the database is strain-level, so the table includes strain. The optional genus and species columns further constrain matching.

Apart from the required columns you can add any other columns to the table that you like. For instance, metadata for the samples, cohort subgroups, notes, etc. Here we have two additional columns tracking the number of metabolites and reactions for each model. If you only want to analyze a single sample the sample_id columns would be the same for each taxon.

The growth medium

A medium must contain reaction and flux columns. Add sample_id when different samples need different media; without it, MICOM applies the same medium to every sample.

from micom.data import test_medium
from micom.qiime_formats import load_qiime_medium

medium = load_qiime_medium(test_medium)
medium

	reaction	flux	metabolite
reaction
EX_glc__D_m	EX_glc__D_m	10.000000	glc__D_m
EX_nh4_m	EX_nh4_m	4.362240	nh4_m
EX_o2_m	EX_o2_m	18.579253	o2_m
EX_pi_m	EX_pi_m	2.942960	pi_m

Building community models

To build one community model for each sample, provide an abundance/taxonomy table and a model database. MICOM databases can be downloaded from the MICOM Zenodo record, or you can provide your own database. A database is optional when every row in the taxonomy table includes a file column pointing to an individual model.

The Batch constructor validates the taxonomy and medium immediately. This makes input errors visible before the more expensive build step.

from micom import Batch
from micom.data import test_medium, test_db

batch = Batch(
    taxonomy=data,
    medium=test_medium,
    model_db=test_db,
)
batch

Batch [mercer_proteus]

4 samples	4 taxa
model DB	✅	strain_models.qza
medium	✅	4 compounds, total flux = 35.88 mmol/(gDW·h)
build	❌
simulated	❌
tradeoffs	❌

cutoff and threads are configured on batch.config rather than passed to Batch.build() directly. The default cutoff keeps taxa representing at least 0.01% of a sample. The build workflow is parallelized using batch.config.threads; set batch.config.build.force_rebuild = True when existing pickle files should be replaced.

The build() workflow returns a manifest and stores it on batch.build_manifest. It also sets batch.out_folder, which later workflows use to locate the community models.

manifest = batch.build("models")
batch

Batch [mercer_proteus]

4 samples	4 taxa
model DB	✅	strain_models.qza
medium	✅	4 compounds, total flux = 35.88 mmol/(gDW·h)
build	✅	found 100.00 ± 0.00% of abundance
simulated	❌
tradeoffs	❌

This will propagate information from your input table as well as give you metrics on how well the samples were matched to the database. Our database only include models for the first 3 E. coli species so you see the workflow could only match 3/4 taxa for each sample. Probably the most important column is the found_abundance_fraction. This one tells which fraction of the sample abundance was matched to the database. You usually want this column to be above 0.5 so the majority of the sample is matched. A value of 1.0 would be perfect but is usually hard to achieve. Values around 0.8 are usually pretty good.

The file column denotes the filename for the built community within the folder specified as out_folder before. You can use the load_pickle function to read individual models and run custom analyses.

from micom import load_pickle

com = load_pickle("models/sample_1.pickle")
print(len(com.reactions))

400

Simulating growth

After build() has completed, grow() simulates each community using batch.medium and the simulation settings in batch.config. The default strategy is cooperative tradeoff. A medium describes which external metabolites can be imported and their maximum fluxes; it can come from a DataFrame, CSV file, or Qiime 2 artifact.

A medium must contain reaction and flux columns. Add sample_id when different samples need different media; without it, MICOM applies the same medium to every sample.

from micom.data import test_medium
from micom.qiime_formats import load_qiime_medium

medium = load_qiime_medium(test_medium)
medium

	reaction	flux	metabolite
reaction
EX_glc__D_m	EX_glc__D_m	10.000000	glc__D_m
EX_nh4_m	EX_nh4_m	4.362240	nh4_m
EX_o2_m	EX_o2_m	18.579253	o2_m
EX_pi_m	EX_pi_m	2.942960	pi_m

The growth medium is represented as a DataFrame with reaction and flux columns. The flux value is the upper bound for import, usually in mmol/gDW/h.

The tradeoff setting controls how much community growth is preserved while individual growth is optimized. It must be between 0 and 1: 1.0 preserves maximum community growth, while smaller values allow more individual optimization. The default value is read from batch.config.simulation.tradeoff.

results = batch.grow()
batch

Batch [mercer_proteus]

4 samples	4 taxa
model DB	✅	strain_models.qza
medium	✅	4 compounds, total flux = 35.88 mmol/(gDW·h)
build	✅	found 100.00 ± 0.00% of abundance
simulated	✅	μᵢ = 0.25 ± 0.20 1/h
tradeoffs	❌

grow() returns a GrowthResults object with three DataFrames: growth_rates, exchanges, and annotations. The same object is stored in batch.results, so it remains available for later analysis and visualization. The has_results property reports whether growth has been simulated.

Choosing a tradeoff parameter

Results depend strongly on the tradeoff parameter. Even though values between 0.3-0.6 usually work well we recommend to run a tradeoff analysis to choose the best parameters for your data set and protocol. If you have already analyzed many samples in your lab and found a particular value to work well in general you may just use that but you should at least run this analysis once. In our paper we found that tradeoff best reproducing in vivo growth rates is the largest tradeoff that allows the majority of the bacteria to grow. Thus, the best tradeoff value is the value providing the best compromise between individual and cooperative growth.

The tradeoff workflow will run growth simulations with several tradeoff values and return the results.

tradeoff = batch.tradeoff()
batch

Batch [mercer_proteus]

4 samples	4 taxa
model DB	✅	strain_models.qza
medium	✅	4 compounds, total flux = 35.88 mmol/(gDW·h)
build	✅	found 100.00 ± 0.00% of abundance
simulated	✅	μᵢ = 0.25 ± 0.20 1/h
tradeoffs	✅	growing fraction: 100.00 - 100.00%

As we see it returns growth rates for each taxon in each sample for various tradeoff values. There is also a tradeoff value of NaN which means optimization of the pure community growth rate without regularization which usually has very bad performance and is provided as a reference. To choose a good value we can count how many of the taxa can grow (growth rate > 1e-6) on average for each of the tradeoff values across all samples.

tradeoff.groupby("tradeoff", dropna=False).apply(
    lambda df: (df.growth_rate > 1e-6).sum()).reset_index()

/tmp/ipykernel_202114/632618865.py:1: FutureWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
  tradeoff.groupby("tradeoff", dropna=False).apply(

	tradeoff	0
0	0.1	16
1	0.2	16
2	0.3	16
3	0.4	16
4	0.5	16
5	0.6	16
6	0.7	16
7	0.8	16
8	0.9	16
9	1.0	16
10	NaN	11

Here all taxa can grow for every tested tradeoff because the example medium provides excess nutrients. For real data, the fraction of growing taxa often declines as the tradeoff increases. The `plot_tradeoff <viz.html>`__ visualization can show this relationship in more detail.

Finding and completing growth media

minimal_medium() finds a medium that enables the requested community and taxon growth rates while minimizing a chosen criterion. It starts from the full set of possible imports.

complete_medium() instead augments an existing medium for the built communities. Use it when you have a plausible starting medium and want MICOM to add only the imports required to meet the growth targets. The two methods account for different scopes: minimal_medium() searches from scratch, while complete_medium() preserves the supplied medium.

import pandas as pd

candidate_medium = pd.DataFrame({"reaction": ["EX_glc__D_m", "EX_o2_m"], "flux": [10, 20]})
candidate_medium

	reaction	flux
0	EX_glc__D_m	10
1	EX_o2_m	20

We can now ask MICOM to complete this medium by adding the smallest amount of overall flux so that all taxa in the database can grow with a growth rate off at least 0.1 1/h.

batch.complete_medium(community_growth=0.1, taxa_growth=0.01, medium=candidate_medium)

	reaction	metabolite	description	flux
0	EX_glc__D_m	glc__D_m	D-Glucose	10.00000
1	EX_nh4_m	nh4_m	Ammonium	0.54528
2	EX_o2_m	o2_m	O2	20.00000
3	EX_pi_m	pi_m	Phosphate	0.36787

So we see that we can achieve growth by adding import for phosphate and ammonium.