Metadata-Version: 2.4
Name: mptool
Version: 1.0.1
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 [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 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).
