Metadata-Version: 2.4
Name: mptool
Version: 1.1.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>=8; 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 [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`) [1]. It can optionally bias the ordering with per-reaction weights, delete in chunks, and favour reactions not yet seen in a sampled MP (`adaptive=True`). Sampled MPs favour MPs with few reactions of their own; see [Sampling](#sampling). 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 CLI prints each MP as it is found, essential reactions included, so the
printed MPs are the ones the run returns. Without `--export` they are the only
output, which is why `--no-verbose` requires `--export`. Their reaction IDs 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.

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.
# Can miss MPs when the fluxes an MP needs are small relative to its
# reactions' upper bounds; it then warns and returns complete=False.
# The iterative method is exact.
mps, mcs, complete = mpt.enumerate_mps(model, subset=subset,
                                       method='direct', graph=False)

# Lossless FVA-based compression is on by default, so the MPs above use
# compressed reaction IDs (see Back-translation below).  Without it:
mps, mcs, complete = mpt.enumerate_mps(
    model, subset=subset, preprocess_kwargs={'compress': False})
```

### Sampling

When complete enumeration is infeasible, random sampling is more practical.
`sample_mps` runs until either `n_samples` distinct MPs have been collected or
`max_t` seconds have elapsed (at least one is required). Duplicates do not
count toward `n_samples`, so if the model has fewer distinct MPs than that —
or fewer that can be drawn: a reaction with weight 0 can make MPs undrawable —
a run without `max_t` never ends. Set `max_t` unless you know there are
enough. Interrupting a run loses the MPs found so far. `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: reactions not yet seen in a sampled MP are kept longer,
# so every reaction that occurs in some MP is found sooner.
mps, acceptance = mpt.sample_mps(model, n_samples=0, adaptive=True, max_t=60)
```

Sampled MPs are not a uniform sample of the MP space. A draw is the MP that
avoids the earliest reactions it can along the ordering, so MPs with few
reactions of their own are drawn more often, and so, for example, are MPs
that share few reactions with others. Frequencies counted per draw are
therefore biased. `mps` holds each distinct MP once, and frequencies over the
distinct MPs found approach the true ones as more of them are found. Neither
weights nor `adaptive=True` removes this bias.

**Weights** are relative: reaction `r` is kept longer than `s` with odds
`w_r : w_s`, so any scale works (raw scores such as BLAST bit scores
included) and multiplying every weight by the same factor changes nothing. A
reaction with weight 0, or missing from `weights`, is tried before every
positive-weight reaction and removed whenever possible; list a reaction to
give it a say. An MP is never drawn while another MP's weight-0 reactions
are a proper subset of its own, and every other MP can be drawn, so 0
excludes wherever an alternative exists while a small positive weight only
disfavours. Each direction of a reversible reaction carries the full weight.
Keys that name no reaction in the subset are ignored, but a `UserWarning` is
raised when none of them does.

With `preprocess=True` (the default), compression does not change the
distribution of sampled MPs, though a given seed gives different MPs with and
without it. Without compression, a chain of fully coupled reactions is
removed, or found to be needed, as soon as the first of them is tried. With
it, each subset member of a merged cluster draws its own key from its weight
(every member counts as 1 without `weights`), and the cluster is placed at
the earliest. This needs the mapping that preprocessing returns, so a model
compressed beforehand and sampled with `preprocess=False` orders each merged
reaction as one reaction, and its sampled MPs are distributed differently.

**`adaptive=True`** multiplies the weight of every reaction that is in none
of the distinct MPs found so far by N, the number found, and leaves the
others at their weight (1 without `weights`; for a merged cluster, every
member's weight). The boost starts with the second distinct MP, and a weight
of 0 stays 0. Once every reaction that occurs in some MP has been seen,
sampling continues as without `adaptive`. `conv_tol` and `conv_window` are
deprecated: their check stops at a roughly fixed number of distinct MPs
rather than when estimates settle.

### Back-translation

The returned MPs and MCSs use preprocessed reaction IDs: with compression
(the default) a flux-coupled cluster is one reaction such as
`'ACONTa+ACONTb'`, and with or without it the reverse half of a reversible
reaction is `'X_rev'`. `return_mapping=True` appends the
preprocessed→original mapping to the return tuple. Each preprocessed reaction
maps to a `{original_id: coefficient}` dict; signs encode direction.

```python
mps, mcs, complete, mapping = mpt.enumerate_mps(
    model, subset=subset, return_mapping=True)

orig_mps = mpt.translate_mps(mps, mapping, subset=subset)
orig_mcs = mpt.translate_mcs(mcs, mapping, subset=subset)
```

Pass the same `subset` you enumerated over. The translation then names
reversed reactions as `X_rev`, as in the preprocessed MPs, and leaves out
reactions outside the subset that compression merged into a cluster with
subset reactions. Expanding `mapping[r]` by hand does neither.

## 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. For sampling, note that `sample_mps` with
`preprocess=False` has no mapping, so on a compressed model it orders each
merged reaction as one reaction, and the sampled MPs are not distributed as
without compression (see [Sampling](#sampling)). Pass `compress=False` to
`preprocess_model` if that matters.

## 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 preprocessed namespace

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:

### Example

```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,
)
```

### 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.**

LP feasibility has this automatically — zeroing bounds only shrinks the feasible
region — but an arbitrary callback need not. Criteria of the form "flux through
R must be at least ...", "the model must still grow", or "this pathway must
remain usable" are safe. Criteria that *reward* inactivity are not: an upper
bound on a flux, such as `EX_ac_e.X <= 5.0`, can fail early and start passing
once more reactions are deleted.

Violate it and returned sets may not be support-minimal, and `chunks=None` and
`chunks=k` may disagree. There is no post-hoc remedy — checking every subset is
exponential.

## Publications

### *Scalable enumeration and sampling of minimal metabolic pathways for organisms and communities*

Cell Systems (2026) — [doi:10.1016/j.cels.2026.101653](https://doi.org/10.1016/j.cels.2026.101653)

Code and notebooks reproducing the results are on the
[`publication/cell-systems-2026`](https://gitlab.com/csb.ethz/mptool/-/tree/publication/cell-systems-2026/paper)
branch, in the `paper` folder.

## 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).
