Metadata-Version: 2.4
Name: mptool
Version: 1.0.0
Summary: Enumeration and sampling of minimal pathways in metabolic (sub)networks.
Home-page: https://gitlab.com/csb.ethz/mptool
Author: Ove Øyås
Author-email: ove.oyas@medisin.uio.no
License: MIT
Project-URL: Source, https://gitlab.com/csb.ethz/mptool
Project-URL: Bug Tracker, https://gitlab.com/csb.ethz/mptool/-/issues
Keywords: metabolism,metabolic networks,metabolic pathways,minimal pathways,pathway analysis,pathway enumeration,pathway sampling,constraint-based modeling,metabolic modeling,systems biology
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Programming Language :: Python :: 3
Classifier: Natural Language :: English
Classifier: Operating System :: OS Independent
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: cobra>=0.29.1
Requires-Dist: gurobipy>=12.0.0
Requires-Dist: networkx>=2.4
Requires-Dist: numpy>=1.21.0
Requires-Dist: scipy>=1.7.3
Provides-Extra: test
Requires-Dist: pytest>=4.6.9; extra == "test"
Provides-Extra: docs
Requires-Dist: sphinx; extra == "docs"
Requires-Dist: furo; extra == "docs"
Provides-Extra: dev
Requires-Dist: mptool[test]; extra == "dev"
Requires-Dist: mptool[docs]; extra == "dev"
Requires-Dist: ruff; extra == "dev"
Requires-Dist: build; extra == "dev"
Requires-Dist: twine; extra == "dev"
Dynamic: license-file

![mptool](https://gitlab.com/csb.ethz/mptool/-/raw/v1.0.0/mptool.png)

# Enumeration and sampling of minimal pathways in metabolic (sub)networks
[![PyPI version](https://badge.fury.io/py/mptool.svg)](https://badge.fury.io/py/mptool) [![Documentation Status](https://readthedocs.org/projects/mptool/badge/?version=latest)](https://mptool.readthedocs.io/en/latest/?badge=latest)

This repository contains a Python package with implementations of methods for metabolic pathway enumeration and sampling as well as notebooks that reproduce the results presented in [our paper](https://doi.org/10.1016/j.cels.2026.101653) [1].

Minimal pathways (MPs) are minimal sets of reactions that need to be active (have non-zero flux) in a metabolic (sub)network to satisfy all constraints on the network as a whole [1]. They can also be defined as the set of *support-minimal* flux patterns from elementary flux vectors (EFVs) [2].

An MP can be found by direct minimization of a mixed-integer linear program (MILP) or by iterative minimization of multiple linear programs (LPs). Both approaches are implemented for enumeration (`enumerate_mps`): the iterative method (default) additionally enumerates minimal cut sets (MCSs), using either a binary integer program (`mcs_method='bip'`, default) [3] or a Berge minimal-hitting-set enumerator (`mcs_method='berge'`). Iterative enumeration can be accelerated by a graph defined by the known MPs to predict unknown MPs (`graph=True`) and by deleting reactions in multiple-reaction chunks rather than one at a time (`chunks`). Random sampling of MPs (`sample_mps`), for when complete enumeration is infeasible, is always iterative; by default it draws uniform random orderings and deletes a single reaction at a time (`weights=None`, `chunks=None`), and can optionally bias the ordering with per-reaction weights, delete in chunks, and update the weights adaptively to reduce length bias [1]. Sampling returns the sampled MPs together with an acceptance rate (the fraction of completed draws — those producing a valid MP — that yield a previously-unseen MP).

## Requirements

To run the code in this repository you will need Python (≥3.9) and the [Gurobi Optimizer](https://www.gurobi.com/) (≥12.0.0). A Gurobi license is required (an academic license is sufficient).

## Installation
You can install `mptool` using [pip](https://pypi.org/project/mptool/):
```
pip install mptool
```
This will also install the minimal requirements `cobra`, `gurobipy`, `networkx`, `numpy`, and `scipy`. The test dependency `pytest` is available via the optional `test` extra (`pip install mptool[test]`).

## Command-line usage

The CLI exposes both entry points as subcommands:

```
python -m mptool enumerate <model file> [--subset FILE] [--bounds FILE] [options]
python -m mptool sample <model file> --n-samples N [--weights FILE] [options]
```

A model in `.xml`, `.json`, or `.mat` format is required. `--subset` accepts a
text file with one reaction ID per line restricting which reactions can appear
in MPs (the candidate set); `--bounds` accepts a CSV file (one
`reaction_id,lb,ub` per line) overriding model bounds for the listed reactions.
The `sample` subcommand additionally accepts `--weights` (a CSV
`reaction_id,weight` per line, in the original namespace) and requires at least
one of `--n-samples` or `--max-t` to terminate.

On `enumerate`, `--method direct` requires `--no-graph` — graph acceleration
applies to the iterative method only — and is not compatible with `--chunks`:

```
python -m mptool enumerate <model file> --method direct --no-graph
```

See `python -m mptool enumerate --help` and `python -m mptool sample --help`
for the full option lists.

### Back-translating compressed CLI output

The reaction IDs the CLI prints are in mptool's *preprocessed* namespace, not
your input model's. With compression on (the default) these include fused
flux-coupled clusters (`A+B`); reversible candidate reactions are also split
into a forward half and a reverse half (`X_rev`), regardless of compression.
Bare stdout output is always in this preprocessed namespace.

To map these IDs back to your original model, pass `--export`: the run bundle
then contains `mapping.csv` — a headerless `new_id,orig_id,coef` table —
alongside `mps.txt`/`mcs.txt`, covering every ID that can appear in the output.
Load it with `mptool.load_mapping`, or translate results directly with
`mptool.translate_mps` / `mptool.translate_mcs`.

## Example usage

The examples below use the e_coli_core model, available from BiGG [4].

```python
import mptool as mpt

# Load model and set a minimal growth rate requirement
model = mpt.load_cobra_model('e_coli_core.xml')
model.reactions.BIOMASS_Ecoli_core_w_GAM.lower_bound = 0.1
```

### Enumeration

All combinations below enumerate the same set of MPs (and MCSs where
applicable) in the model's boundary reactions — every minimal combination of
metabolite uptakes and secretions that supports growth. This subset is small
enough that enumeration finishes within seconds.

```python
subset = model.boundary

# Iterative LP-based enumeration (default) with graph-based prediction
# and BIP minimal-cut enumeration.
mps, mcs, complete = mpt.enumerate_mps(model, subset=subset)

# Iterative without graph acceleration.
mps, mcs, complete = mpt.enumerate_mps(model, subset=subset, graph=False)

# Iterative with the Berge minimal-hitting-set enumerator for MCSs
# (alternative to the default 'bip').
mps, mcs, complete = mpt.enumerate_mps(model, subset=subset,
                                       mcs_method='berge')

# Direct MILP minimization: finds one MP per solve; no MCSs returned.
mps, mcs, complete = mpt.enumerate_mps(model, subset=subset,
                                       method='direct', graph=False)

# Lossless FVA-based compression applied before enumeration.  Returned
# MPs use compressed reaction IDs; back-translate via the mapping.
mps, mcs, complete, mapping = mpt.enumerate_mps(
    model, subset=subset,
    preprocess_kwargs={'compress': True, 'tighten': True},
    return_mapping=True,
)
```

### Sampling

When complete enumeration is infeasible, random sampling is more practical.
`sample_mps` runs until either `n_samples` MPs have been collected or `max_t`
seconds have elapsed (at least one is required). `max_t` budgets the sampling
loop only — preprocessing runs before the clock starts, so total runtime is
preprocessing plus `max_t`. The same holds for `max_t` in `enumerate_mps`. It
returns a 2-tuple
`(mps, acceptance)`: `mps` is a list of MPs in sampling order (so reaction
frequencies and convergence statistics can be recomputed directly from it)
and `acceptance` is the acceptance rate — the fraction of completed draws
that yielded a new MP rather than a duplicate (a low rate signals the MP
space is nearly exhausted; prefer `enumerate_mps` then). The full model is
used here (no `subset` argument).

```python
import random

# Uniform sampler (default): single-reaction-deletion ordering, no weights.
mps, acceptance = mpt.sample_mps(model, n_samples=100)

# Chunked deletion (faster per sample on large subsets; degrades gracefully
# to single-reaction deletion when fewer than `chunks` candidates remain).
mps, acceptance = mpt.sample_mps(model, n_samples=100, chunks=4)

# Weighted sampler: prior weights bias the random ordering
# (higher weight → reaction more likely to be retained in sampled MPs).
weights = {r.id: random.random() for r in model.reactions}
mps, acceptance = mpt.sample_mps(model, n_samples=100, weights=weights)

# Adaptive sampler: weights are updated after each sample to reduce
# length bias and improve coverage of the MP space.
mps, acceptance = mpt.sample_mps(model, n_samples=0, adaptive=True, max_t=60)
```

### Back-translation when `compress=True`

With `compress=True` the returned MPs use compressed reaction IDs (e.g.
`'ACONTa+ACONTb'` for a flux-coupled cluster, or `'X_rev'` for the reverse
half of a reversible reaction). `return_mapping=True` appends the
preprocessed→original mapping to the return tuple. Each merged reaction maps
to a `{original_id: coefficient}` dict; signs encode direction.

```python
mps, mcs, complete, mapping = mpt.enumerate_mps(
    model, subset=subset,
    preprocess_kwargs={'compress': True, 'tighten': True},
    return_mapping=True,
)

# Expand a single MP back to original reaction IDs (union over members):
original_ids = {orig for r in next(iter(mps)) for orig in mapping[r]}
```

## Preprocessing

`enumerate_mps` and `sample_mps` preprocess the model internally — copying it,
applying user bounds, removing flux-inconsistent reactions (fastercc),
optionally applying FVA-based bound tightening and lossless compression, and
making the subset forward-irreversible. The same pipeline is exposed as a
standalone function:

```python
preprocessed, subset, mapping = mpt.preprocess_model(
    model, subset=['BIOMASS_Ecoli_core_w_GAM', 'EX_glc__D_e'],
    bounds={'BIOMASS_Ecoli_core_w_GAM': (0.1, 1000.0)},
    compress=True, tighten=True,
)
mps, mcs, complete = mpt.enumerate_mps(
    preprocessed, subset=subset, preprocess=False,
)
```

This is useful when the same model is reused across multiple
enumeration / sampling calls — preprocessing is the expensive step and
running it once amortises its cost.

## Writing requirements functions

The `requirements` parameter of `enumerate_mps` / `sample_mps` accepts a
callable that imposes additional feasibility criteria on top of LP
feasibility. It is a **top-level keyword argument** and requires
`preprocess=False`: preprocess the model yourself via `preprocess_model` and
pass the result in. (It is not compatible with `method='direct'`.)

The callable has signature `(model) -> (bool, set[str])`; a plain `bool`
return is also accepted (treated as `(bool, set())`). On each invocation it
receives a fresh independent copy of the current Gurobi model and may freely
mutate that copy without affecting the search. The copy is not solved before
the callback runs, so a callback that inspects solution values must call
`model.optimize()` itself. It returns a 2-tuple `(feasible, zero_flux_ids)`:

- `feasible` (`bool`): True if the active reactions satisfy the criterion.
- `zero_flux_ids` (`set[str]`): Reaction ID strings known to carry zero flux
  under the current state; these can be pruned without re-solving the LP.
  Return an empty set if no such information is available.

### The monotonicity precondition

Your callback must satisfy one property, and mptool cannot check it for you:

> **Once a set of deactivated reactions makes your criterion fail, deactivating
> more reactions must not make it pass again.**

Equivalently: if your criterion holds for some set of deleted reactions, it must
also hold for every smaller set of deleted reactions. LP feasibility has this
property automatically — forcing a reaction's bounds to zero only shrinks the
feasible region, so feasibility can never be *recovered* by deleting more. An
arbitrary callback need not.

The search relies on it in two places:

- **Support-minimality.** Reactions are deleted greedily along an ordering. When
  deleting reaction `x` breaks the criterion, `x` is reinstated and kept in the
  MP — but that decision is made *at that point*, with only the reactions before
  `x` in the ordering deleted. More get deleted afterwards. Concluding that `x`
  is still needed at the end is exactly the property above. Without it, the
  returned set can contain reactions that were removable by the time the search
  finished.
- **`chunks` equivalence.** Chunked deletion (`chunks=k`) tries to delete a whole
  block at once and only splits it when that fails. Accepting a whole-block
  deletion assumes each individual member was deletable too — the property
  again. Without it, `chunks=None` and `chunks=k` can return different sets from
  the same ordering.

**A callback that violates it.** Anything that *rewards* inactivity — an upper
bound on a flux, or a requirement that some reaction be off:

```python
def bad_requirements(m):           # DON'T: not monotone
    m.optimize()
    if m.status != GRB.OPTIMAL:
        return False, set()
    # "byproduct secretion must stay below 5" — deleting MORE reactions
    # lowers this flux, so a failing state can start passing again.
    return m.getVarByName('EX_ac_e').X <= 5.0, set()
```

Deleting `x` early may leave acetate secretion at 8, so the criterion fails and
`x` is kept. Deleting `y` later drops secretion to 3 — and now `x` was never
needed. The returned set includes `x` and is not minimal.

Criteria of the form "flux through R must be at least ...", "the model must
still grow", "this pathway must remain usable" are monotone and safe: deleting
more reactions can only make them harder to satisfy.

**If your callback cannot be made monotone**, use `validate_mps`, which tests
support-minimality *directly* rather than inferring it from the search order: it
takes each returned MP, forces every subset reaction outside the MP off, and
then forces each MP reaction off one at a time, re-checking your callback each
time. It accepts `requirements` itself, so it gives an honest per-MP verdict
under any callback:

```python
verdicts = mpt.validate_mps(preprocessed, mps, subset,
                            requirements=my_requirements)
bad = [mp for mp, ok in zip(mps, verdicts) if not ok]
```

One caveat on that escape hatch: it checks that no *single* reaction can be
dropped. For a monotone callback that is equivalent to full minimality; for a
non-monotone one it is the weaker, one-at-a-time statement, since dropping two
reactions together could still succeed where dropping either alone fails.

Because `requirements` runs with `preprocess=False`, the callback sees the
namespace of the model you passed in — the output of `preprocess_model`.
That model is forward-irreversible, so every reaction has non-negative flux
bounds: reversible reactions are split into a forward half (original name)
and a reverse half (`_rev` suffix); a reaction `PFK` with `lb=-10, ub=0`
becomes `PFK_rev` with `lb=0, ub=10`. With `compress=True`, IDs may differ
further from the original model (merged reactions use combined `'A+B'` names;
blocked reactions are absent).

Requirements functions must therefore be written against the preprocessed
model. The cleanest approach is to preprocess once, inspect the resulting
variable names, and write the function using those names:

```python
import mptool as mpt
from gurobipy import GRB

model = mpt.load_cobra_model('e_coli_core.xml')
model.reactions.BIOMASS_Ecoli_core_w_GAM.lower_bound = 0.1

# Preprocess once, then inspect the preprocessed reaction IDs.
preprocessed, subset, _ = mpt.preprocess_model(
    model, subset=model.boundary, compress=False)
preprocessed_ids = {v.varName for v in preprocessed.getVars()}

# Write a requirements function using the preprocessed namespace.  The copy
# is not solved for you, so optimize it before reading solution values.
def my_requirements(m):
    m.optimize()
    if m.status != GRB.OPTIMAL:
        return False, set()
    growth = m.getVarByName('BIOMASS_Ecoli_core_w_GAM')
    return (growth is not None and growth.X >= 0.1 - 1e-6), set()

# Pass it as a top-level kwarg alongside preprocess=False.
mps, mcs, complete = mpt.enumerate_mps(
    preprocessed, subset=subset, preprocess=False,
    requirements=my_requirements,
)
```

## Reproducing results from publication

In the `paper` folder we provide code for reproducing the results presented in [our paper](https://doi.org/10.1016/j.cels.2026.101653) [1]. Models were obtained from [BiGG](http://bigg.ucsd.edu) [4] and [Virtual Metabolic Human](https://www.vmh.life/) [5].

### Example network analysis

The notebook `example.ipynb` compares pathways using an example network. The MATLAB script `efm_efv_example.m` uses CellNetAnalyzer [6] to enumerate elementary flux modes and vectors for the example network.

### Benchmarking methods

The notebook `benchmarking.ipynb` analyzes results from benchmarking of methods. The script `benchmarking.py` was used to perform the benchmarking on a cluster.

### *E. coli* core metabolism analysis

The notebook `e_coli_analysis.ipynb` samples and analyzes MPs from *E. coli* core metabolism in the context of the full genome-scale network.

### Host-microbe interaction analysis

The notebook `host_microbe_analysis.ipynb` samples, enumerates, and analyzes MPs in a sequentially constrained host-microbe model of the human gut.

### Butyrate-producing community analysis

The code in the archive `butyrate_analysis.zip` enumerates and analyzes minimal butyrate-producing microbial communities. This also requires the code from the experimental study that was used for comparison to predictions, which can be found [here](https://github.com/RyanLincolnClark/DesignSyntheticGutMicrobiomeAssemblyFunction).

## Citation

If you use `mptool` for a scientific publication please cite [our paper](https://doi.org/10.1016/j.cels.2026.101653) [1].

## Development

`mptool` is developed and maintained by Ove Øyås, with development done in
collaboration with Claude Code.

## References

[1] O. Øyås, A. Theorell, and J. Stelling. "Scalable enumeration and sampling of minimal metabolic pathways for organisms and communities". *Cell Systems* (2026). https://doi.org/10.1016/j.cels.2026.101653

[2] S. Klamt et al. "From elementary flux modes to elementary flux vectors: Metabolic pathway analysis with arbitrary linear flux constraints". *PLoS Computational Biology* 13.4 (2017).

[3] H.S. Song et al. "Sequential computation of elementary modes and minimal cut sets in genome-scale metabolic networks using alternate integer linear programming". *Bioinformatics* 33.15 (2017).

[4] Z.A. King et al. "BiGG Models: A platform for integrating, standardizing, and sharing genome-scale models" *Nucleic Acids Research* 44.D1 (2016).

[5] A. Noronha et al. "The Virtual Metabolic Human database: integrating human and gut microbiome metabolism with nutrition and disease" *Nucleic Acids Research* 47.D1 (2018).

[6] A. von Kamp et al. "Use of CellNetAnalyzer in biotechnology and metabolic engineering" *Journal of Biotechnology*, 261 (2017).
