Metadata-Version: 2.4
Name: rafkit
Version: 0.6.0
Summary: Autocatalytic (RAF) sets in catalytic reaction networks: maximal RAFs, irreducible cores, and Kauffman binary polymer models.
Author: James P. Galasyn, Claude Théodore
License: MIT
Project-URL: Homepage, https://github.com/jimgalasyn/rafkit
Project-URL: Issues, https://github.com/jimgalasyn/rafkit/issues
Keywords: autocatalytic-sets,RAF,origin-of-life,chemical-reaction-networks,systems-chemistry,binary-polymer-model,catalysis,abiogenesis
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Artificial Life
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Provides-Extra: cac
Requires-Dist: scipy>=1.10; extra == "cac"
Provides-Extra: test
Requires-Dist: pytest>=8; extra == "test"
Requires-Dist: pytest-cov>=4; extra == "test"
Requires-Dist: pytest-xdist>=3; extra == "test"
Requires-Dist: scipy>=1.10; extra == "test"
Dynamic: license-file

# rafkit

[![CI](https://github.com/JimGalasyn/rafkit/actions/workflows/ci.yml/badge.svg)](https://github.com/JimGalasyn/rafkit/actions/workflows/ci.yml)
[![codecov](https://codecov.io/gh/JimGalasyn/rafkit/branch/main/graph/badge.svg)](https://codecov.io/gh/JimGalasyn/rafkit)
[![PyPI](https://img.shields.io/pypi/v/rafkit.svg?cacheSeconds=3600)](https://pypi.org/project/rafkit/)
[![Python](https://img.shields.io/pypi/pyversions/rafkit.svg?cacheSeconds=3600)](https://pypi.org/project/rafkit/)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE)
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.21954795.svg)](https://doi.org/10.5281/zenodo.21954795)

Autocatalytic (RAF) sets in catalytic reaction networks — maximal RAFs, irreducible
cores, Kauffman binary polymer models, and interoperability with
[CatReNet](https://github.com/husonlab/catrenet).

Pure Python and NumPy. No Java, no GUI, no install beyond `pip`.

```python
from rafkit import binary_polymer, max_raf, sample_irrraf
import numpy as np

net = binary_polymer(max_len=8, food_len=2, p=1.5e-3, cleavage=True)
raf = max_raf(net)
print(raf.size, "reactions in the maximal RAF")

core = sample_irrraf(net, raf.reactions, np.random.default_rng(0))
print(len(core), "reactions in one irreducible core")
```

## Why this exists

RAF theory (Hordijk & Steel 2004) formalises collectively autocatalytic sets: a set of
reactions is a RAF over a food set when every reaction is catalysed by something the
set can make, and every reactant can be built up from food using the set. The
reference implementation, **CatReNet**, is an excellent Java/JavaFX desktop
application. This is a small library for people who want the same algorithms inside a
Python analysis pipeline.

## Validated against the reference implementation

`tests/data/catrenet_polymer_n6.crs` was generated by CatReNet's own `polymer-tool`,
and the expected counts in `tests/test_crs.py` are what CatReNet's `catrenet-tool`
reports on it. The test suite therefore checks this implementation against an
independent one on every run, with no Java required.

| algorithm | rafkit | CatReNet 1.1.0 |
|---|---|---|
| `max_raf` | 183 | 183 |
| `catrenet_strictly_autocatalytic` | 175 | 175 |
| `max_raf_strict` | 161 | *(different object — see below)* |

**A documented divergence.** CatReNet's `strictlyAutocatalyticMaxRaf` *filters* the
maximal RAF for reactions having a non-food catalyst, without re-refining, so its
result need not itself be a RAF. `max_raf_strict` imposes the same condition inside
the fixpoint, so its result is a RAF, and is correspondingly smaller. Both are
available; they answer different questions. CatReNet's behaviour was reproduced by
black-box inference from its output — no CatReNet source was read or used.

## Calibration

The RAF phase transition in Kauffman's binary polymer model, measured here against
the published value of *f* ≈ 1.20 (Steel, Hordijk & Smith 2012, n=10, t=2), where
*f* = p|R| is the mean number of catalysed reactions per molecule:

| model | transition |
|---|---|
| ligation only | *f* ≈ 4.7 |
| ligation + cleavage | *f* ≈ 3.1 |
| ligation + cleavage, **catalysis paired per reversible reaction** | 0 seeds at *f* ≤ 1.22, all seeds by *f* ≈ 1.59 |

Two conventions have to match before any comparison to the literature means anything:
the model must include **cleavage**, and a reversible cleavage–ligation pair must be
counted as **one** catalysed reaction, not two. Use `net.catalysis_level` — not
`mean_catalysed_per_molecule` — whenever a number is placed beside a published *f*.

## What's implemented

| | |
|---|---|
| `max_raf` | maximal RAF, by fixpoint (Hordijk & Steel 2004) |
| `max_raf_strict` | maximal RAF whose catalysts must be non-food products |
| `catrenet_strictly_autocatalytic` | CatReNet's similarly-named filter, for interop |
| `sample_irrraf` | one irreducible RAF, by randomised shrinking (Steel, Hordijk & Smith 2012) |
| `irrraf_census` | how many *distinct* irreducible cores a network carries |
| `exploitability` | share of RAF products contributing no catalysis back |
| `is_food_catalysed` | whether a core runs on food catalysis alone, and so carries no heredity |
| `core_raf` / `has_unique_irraf` | Huson, Xavier & Steel's polynomial test for a *unique* irreducible RAF |
| `catalytically_reachable` | what can be made without any spontaneous reaction |
| `binary_polymer` | Kauffman binary polymer generator, with optional cleavage |
| `ReactionNetwork` | arbitrary catalytic reaction systems, same protocol |
| `read_crs` / `write_crs` | CatReNet's CRS interchange format |
| `to_pnml` / `write_pnml` | PNML export (ISO/IEC 15909-2) for the Petri net ecosystem |
| `simulate` | Gillespie direct method — watch subRAFs seed themselves into existence |
| `max_urafs` | uninhibited RAFs, when a molecule can prevent a reaction |
| `run_serial_dilution` / `run_cstr` | dilution protocols for growing–dividing compartments — **not a RAF algorithm**, see below |
| `permeation_flux` / `permeable_by_length` | size-selective transport across a compartment membrane — **not a RAF algorithm**, see below |
| `is_thermodynamically_consistent` | can these reactions all run forward at once? — **not a RAF algorithm**, see below |
| `BondEnergies` / `rate_constants` | free energy of a polymer chemistry, and the rate constants it forces — **not a RAF algorithm**, see below |
| `detailed_balance_residual` | whether a set of rate constants is consistent with the free energies |
| `elongation_ratio` / `mean_length` / `sequence_correlation_length` | the equilibrium ensemble, closed form |
| `Kinetics` / `kinetics_from_energies` | rate constants a simulator can run on — `simulate(..., kinetics=...)` |
| `unpaired_catalysis` | reversible pairs whose two directions have different catalysts, which is impossible |

Every algorithm carries hand-computed known-answer tests, because a RAF algorithm that
is subtly wrong produces plausible numbers rather than errors.

## Four modules are deliberately off-theme: `dilution`, `permeation`, `thermo` and `autocatalysis`

Everything above takes a `ReactionNetwork` and asks a RAF question of it. `rafkit.dilution`
takes no network at all — it is a **two-species ordinary differential equation with no RAF
structure**, reproducing the minimal model of Matsubara, Ameta, Thutupalli, Nghe & Krishna
([arXiv:2211.03155](https://arxiv.org/abs/2211.03155)).

It earns its place for one reason: **it is this library's only *analytic* calibration.**
Every other check here is against a reference implementation (CatReNet) or a published
figure (Steel, Hordijk & Smith) — matched to a count, or to a shape. Matsubara et al. derive
closed-form conditions, which can be hit or missed to nine significant figures:

| their claim | status |
|---|---|
| `r(x)x` linear ⇒ only the symmetric trajectory is stable, **no bistability** | reproduced |
| `r(x)x = ε + κx²` ⇒ **bistability** at their `Δt`=1, `κ`=8, `ε`=0.5, `φ`=1 | reproduced |
| bistability **lost above a critical cycle interval** | reproduced |
| their equation (2), **parameter-free** | reproduced to 2×10⁻⁹ (asserted at 1e-8) |

The two parameterisations also agree on the *sign*: their equation (3) makes the
amplification factor crossing 1 the sufficient condition for bistability, and it measures
0.900 for the linear flux against 1.243 for the quadratic — the numerical test and the
analytic criterion picking out the same case from independent computations.

The wider use is that RAF work increasingly runs networks inside growing, dividing
compartments, where the dilution protocol is a modelling choice that changes the answer.
This gives that choice a validated implementation and a benchmark, independent of any RAF
structure. If you only want RAF algorithms, ignore this module; nothing else imports it.

### `permeation` — one line of transport physics, and a trap worth a module

`rafkit.permeation` is smaller and earns its place differently. RAF work increasingly runs
networks inside compartments embedded in a shared medium, and every such model needs a rule
for what crosses the boundary. The rule implemented is Hordijk, Naylor, Krasnogor &
Fellermann's ([*Life* **8**(3), 33, 2018](https://doi.org/10.3390/life8030033)) — *"molecules
are allowed to permeate compartment membranes if their lengths do not exceed a certain
threshold. Permeation is proportional to the concentration difference."*

**"Proportional to the concentration difference" is not "proportional to the count
difference",** and the two coincide only when compartment and medium have the same volume. In
a spatial model they generally do not: in the paper above a compartment of radius 0.5 sits in
a diffusion voxel of 2.5 × 2.5 — and since that world is two-dimensional ("if the world type
is set to 2D then Y is forced to 1"), the voxel is a slab of unit thickness. **Both are then
volumes**: a sphere of 0.524 against a slab of 6.25, a ratio of **11.9**. Quoting a voxel
*area* against a sphere *volume* would be dimensionally meaningless. Writing the flux as
`P · (n_out − n_in)` silently asserts they are the same size, and produces plausible numbers
rather than an error — which is the failure mode this whole library is written against.

How much it matters: reproducing that paper's own induction experiment with every printed
parameter taken from [the authors' published input files](http://ico2s.org/data/extras/compartments/),
the count-difference form **cannot match both published arms at any permeability** — swept, it
reaches the control value at an effect ratio of 1.15 on one branch or 2.44 on the other,
bracketing the published 1.60 without hitting it. The concentration form reproduces both arms
(16.5 ± 5.9 against their 16.3; 27.8 ± 6.0 against their 26.0) with a single free parameter.

⚠ Note the calibration tier: that is a published *figure* matched with one fitted parameter, so
it sits with this library's reference-implementation and figure checks — **not** with
`dilution`, which remains the only analytic anchor here. The module's own tests are
deterministic properties (equal concentrations give zero flux; equal *counts* do not; flux
bounded by what is present), because a library gate should be fast and exact; the stochastic
reproduction lives in the downstream research client.

### `thermo` — the free energy the rest of the library leaves implicit

`rafkit.thermo` is the third off-theme module, and it is here because **the rest of this
library already has a thermodynamics — an unwritten one, and it is the wrong one.**

`gillespie` gives every reaction a unit rate constant. In a cleavage–ligation chemistry that
means `k_f = k_r` for every reaction, so `K_eq = 1` and `ΔG° = 0`: **every polymer is
isoenergetic with the parts it is made of.** Nothing is more stable than anything else, no
sequence is preferred over any other, and the polymer growth those runs show is driven
entirely by the food boundary condition. That is not a badly chosen parameter; it is free
energy being *absent* while the model talks as though it were present. The check is one call —
against a bond energy of −1 to −3 with an association cost of 0.5, unit rate constants score a
`detailed_balance_residual` of 2.5, and the rates this module builds score 2×10⁻¹⁶.

**Catalysis has the same problem one level down.** Where the uncatalysed rate is zero,
"catalysis" is not acceleration but *enablement*: the catalyst decides whether the reaction
exists at all, so "which species catalyses what" becomes a choice of which reactions there
are, wearing a kinetic name. A catalyst that enables also moves the equilibrium — from
unreachable to reachable — and no catalyst does that. Here a catalyst is a ratio applied to
**both directions at once**, which is the only form leaving `K_eq` alone; the residual for
enablement is `inf`, not a large number. For scale, 100× is a barrier drop of 11.4 kJ/mol at
298 K, about one hydrogen bond.

Detailed balance is structural rather than imposed: the barrier is split between the two
directions in the Brønsted way, `ΔG‡_f = barrier + β·ΔG` and `ΔG‡_r = barrier − (1−β)·ΔG`, so
`k_f/k_r = exp(−ΔG/RT)` identically for **every** barrier and every `β`. There is no parameter
setting that violates it, and `k_uncat` stops being a switch: fix the barrier and both rates
follow, neither of them zero.

**Why three bond energies and not one.** A uniform bond energy gives every equal-length
sequence identical free energy, so no thermodynamic sequence preference can exist — in a unary
alphabet the question cannot arise, in a binary one it is the whole point. But "uniform" is not
the tight condition. The ensemble these energies induce is exactly a **one-dimensional Ising
chain**, whose transfer matrix loses its second eigenvalue whenever

```
ε = E₀₀ + E₁₁ − E₀₁ − E₁₀ = 0
```

— the *additive* case `E[a][b] = h(a) + g(b)`, where the bond energy says something about the
left residue and something about the right one and nothing about the pair. (Under
`BondEnergies.symmetric`, where `E₀₁ = E₁₀`, that is the familiar `E₀₀ + E₁₁ − 2·E₀₁`. The
matrix is not required to be symmetric — a directional backbone need not be — and the doubled
form reports a preference that is not there when it is not.) Every additive assignment has a sequence correlation length
of exactly zero, uniform or not, so three energies satisfying `E₀₁ = (E₀₀+E₁₁)/2` buy nothing
over one. **It is the non-additivity `ε` alone that makes ordering thermodynamically visible**:
`ε > 0` favours alternation, `ε < 0` favours blocks. That is where a thermodynamic basis for
templating would have to come from, as opposed to an imposed rule.

Two exact anchors, both in `tests/test_thermo.py`:

| claim | status |
|---|---|
| ΔG of a ligation is the **junction bond alone**, and every split of every sequence agrees | exact — this is Wegscheider consistency |
| equilibrium length distribution is geometric with ratio `ρ = λ_max(T)`, mean `1/(1−ρ)` | exact for an additive assignment, from the first **bond** |

⚠ That last row was wrong in its first form, and the correction is worth stating: the geometric
law starts at the first *bond*, not the first *molecule*. A monomer has no bonds, so the step
from length 1 to length 2 is a boundary term and need not equal `ρ`. A uniform assignment hides
this — its first step happens to equal `ρ` — so the null case is the one case where the effect
is invisible. On an additive-but-not-uniform assignment at `ρ = 0.4` the first step is 0.377 and
`mean_length` overstates the true number-average by 1.4%.

⚠ And a third, found only because the binary case is the easy one: **Perron–Frobenius
constrains the leading eigenvalue and nothing else.** A positive matrix of size 3 or more may
have complex subdominant eigenvalues — a randomly drawn 4×4 bond-energy matrix does on the
first try — so rejecting them as impossible refused an ordinary nucleotide chemistry. They are
legitimate: a complex pair is a correlation that *oscillates* as it decays, and the decay
length is set by the modulus either way.

**Reaching the simulator.** `kinetics_from_energies` turns bond energies into a `Kinetics` —
per-reaction *uncatalysed* rate constants plus the factor a **present** catalyst applies — and
`propensities(..., kinetics=...)` runs on it. The split is deliberate: `thermo` says what the rate
constants are, and whether a catalyst is present at any instant is state, so the simulator never
holds a "catalysed rate constant" it could apply to one direction.

⚠ **The existing model was already in this form and did not say so.** `gillespie`'s uncatalysed
factor of 20 *is* `Kinetics.uniform(net.n_reactions, 1/20, 20)`, and reproduces its propensities to the
last ulp — so catalysis was never the missing piece. What was missing is the equilibrium: unit
constants make `k_f = k_r` whether or not a catalyst is present, so `K_eq = 1` regardless.

The consequence that can change a trajectory: a chemistry built this way has its **stationary point
at the thermodynamic equilibrium.** Balance is on the *combinatorial factors* —
`k_f · combos_forward = k_r · combos_reverse` — which for `a + b → ab` with `a ≠ b` is the familiar
`n_ab/(n_a·n_b) = K`. A present catalyst raises both directions by the enhancement **without moving
that point**. Under unit constants the balance sits at `K = 1` for every reaction, whatever the
molecules are.

⚠ **Not so for a self-ligation.** `a + a → aa` takes the pair count `n_a(n_a−1)/2`, so it balances at
`n_aa = K·n_a(n_a−1)/2` and the count ratio `n_aa/n_a²` tends to **`K/2`**. Measured on `0 + 0 → 00`
under unit constants at `n_a = 20`, balance is at `n_aa = 190`, not 400. The factor is the standard
stochastic symmetry number and `_pair_count` is correct — but **the map from ΔG to a count ratio is
not uniform across the chemistry**, and anyone reading equilibrium constants off a trajectory needs
the qualification.

⚠ `binary_polymer(paired_catalysis=False)` is not a variant chemistry — it is a **thermodynamically
impossible** one. Drawing the two directions' catalysts separately gives molecules that accelerate a
ligation but not its cleavage, and whenever such a molecule is present the reaction is a free-energy
source; measured on a hand-built case, 100× net flux from nothing. `unpaired_catalysis` finds them
and `kinetics_from_energies` refuses them.

⚠ A `simulate` run with `kinetics` is still **driven**, not closed: the food floor holds a chemical
potential at the boundary. It does not relax to the equilibrium ensemble computed above, and
comparing the two directly would be comparing a driven steady state to an equilibrium.

**`dg_assoc` is required, with no default.** The association cost — the standard-state price of
turning two molecules into one — is the only sequence- and length-independent term in a ligation,
and `0.0` is not a neutral absence but the claim *"joining is free"*. For scale, Ross & Deamer
(*Life* **6**(3):28) put phosphodiester formation at **+3.3 kcal/mol at 85 °C (≈ +4.7 RT,
K₁ ≈ 1e-3)**, so zero is not a small value of this quantity — it is a different claim. It also
decides whether an equilibrium exists at all: at monomer 0.5 with `E = −1`, `dg_assoc = 0` puts the
elongation ratio **above 1** (runaway, no equilibrium) and `+4.7` brings it back below.

Same rule as `permeation_flux` requiring both volumes: where a value that looks like an absence is
really an assertion, there is no default that would let it happen silently. ⚠ The single exemption
is `sequence_correlation_length`, where the argument **provably cancels** — requiring a value that
cannot change the answer only trains the reflex that makes the requirement worthless everywhere
else.

⚠ And it is where **water activity** enters, since a ligation releases water: mass action for
`N_m + N_n ⇌ N_{m+n} + H₂O` puts `+RT·ln(a_W)` in exactly this term, negative when dry. So the old
default silently asserted `a_W = 1` — permanently wet. **Not implemented**: that is the form, not a
number.

⚠ Calibration tier: **algebraic**, not empirical. Everything above is an identity that holds or
does not, so it sits beside `dilution` rather than beside the figure reproductions — but it
reproduces no experiment and calibrates against no published number. It says the model is
*consistent*, not that it is *right*.

### `autocatalysis` — the one place this library can be told it is wrong by someone else

`rafkit.autocatalysis` decides a question RAF theory cannot pose: **given that these reactions must
run in these directions, does any assignment of chemical potentials make that happen?** That is the
CAC question of Kosc, Kuperberg, Rajon & Charlat, [*PNAS* **122**(18) e2421274122
(2025)](https://doi.org/10.1073/pnas.2421274122).

**The whole module is one reduction.** With `x = e^μ` and barrier factor `b_i = e^{−G‡_i}`, local
detailed balance gives `v_i = b_i(∏x^{S⁻} − ∏x^{S⁺})`. Since `b_i > 0`, **the barrier scales the flow
but cannot flip its sign**, so with `y = ln x = μ` reaction `i` runs forward exactly when
`(Sᵀy)_i < 0`. The question becomes **strict linear feasibility of `Sᵀ y < 0`** — a linear program in
the chemical potentials, and nothing else.

⚠⚠ **The verdict therefore depends on no rate constant and no barrier**, which the paper states
outright. The API accepts none, and a test asserts that it accepts none.

**Two independent methods, required to agree.** Gordan's theorem says exactly one of `Sᵀy < 0` (a
witness) and `S w = 0, w ≥ 0, w ≠ 0` (a certificate) can hold. Both are computed and disagreement
**raises** rather than picking a winner.

**Checked against the published answers** — Kosc's Fig. 4, two cores sharing `{R2, R3}`:

| | verdict | source |
|---|---|---|
| `{R1,R2,R3,R4}` alone | consistent | **Theorem 2** — a single PAC always is |
| `{R1',R2,R3,R4'}` alone | consistent | **Theorem 2** |
| both together | **inconsistent** | **Box 2** — a multiPAC that is not a multiCAC |

⚠ That network is *reconstructed* from Box 2's flow equations (the paper draws it as a figure) and
cross-checked against the figure's composition glyphs: every reaction mass-balances, both cores are
autocatalytic in `e4`, and they share exactly two reactions as the caption says.

**And the certificate explains itself.** It comes back with every weight equal — the six reactions at
unit flux return the system to its starting composition. Then `Σ w_i A_i = −(S w)·y = 0`, so the
affinities cannot all be positive: **a cycle that returns to its starting composition cannot be
downhill all the way round.** The second law, as a linear-algebra identity.

**Why it matters more than a bigger test suite.** Theorem 2 is an *external* guarantee — a single PAC
is always consistent — so **every network this library can generate is a pass/fail case adjudicated
by someone else.** It is the one place rafkit can be shown wrong without anyone here noticing first.

⚠ Needs `scipy` for the LP: `pip install rafkit[cac]`. Imported lazily, so `import rafkit` and every
other module stay numpy-only.

## A reaction network is a Petri net

Species are places, reactions are transitions, molecule counts are tokens. `write_pnml`
exports to PNML (ISO/IEC 15909-2), so these networks open in Petri net editors, model
checkers, and the unfolding tools that compute the causal structure of a run.

Three things need care, and each is explicit rather than silent:

- **Catalysis becomes a self-loop.** P/T nets have no read arc, so a catalyst is a pair
  of arcs, consuming the token and returning it.
- **Alternative catalyst sets become separate transitions**, named `r1`, `r1#2`, …, since
  a transition's preset is a conjunction and cannot express "either set".
- **Food gets source transitions**, because RAF food is inexhaustible and no initial
  marking expresses that — a marking of *n* deadlocks after *n* uses.

Inhibition has no `ptnet` representation and is written as a `toolspecific` annotation,
with a warning in the file: a reader that ignores it gets a *different system*.
Reactions requiring a catalyst that nothing provides are omitted and counted, since
emitting them unconstrained would make them freely fireable — the opposite of the intent.

## Catalysis is a relation, not a list

`catalysts[r]` is a set of **alternative catalyst sets**, following Huson, Xavier &
Steel (2024). Any one set being fully present suffices, and each set is a conjunctive
requirement:

| `catalysts[r]` | meaning |
|---|---|
| `{{a}, {b}}` | *a* **or** *b* — the simple case, and what a flat list of catalysts meant |
| `{{a, d}, {e}}` | (*a* **and** *d*) **or** *e* |
| `{}` | **must** be catalysed, and nothing does: never in a RAF |
| `{frozenset()}` | **may proceed uncatalysed**; always satisfied |

The last two rows are a real distinction rather than a technicality — in the §2.4 system
of that paper it decides which reactions can join an RAF — and a flat list collapses
them. In CRS, a braced group is conjunctive: `[{a,d}, e]`, with `[]` and `[{}]` for the
last two rows.

Constructors still accept a plain iterable of molecules and normalise it, so simple
systems stay simple to write.

## Inhibition

A molecule can prevent a reaction. `max_urafs` returns the **uninhibited RAFs** of
Hordijk & Steel (2012), and returns a *collection* rather than one set, because
inhibition destroys the monotonicity that makes a maximal RAF unique — there is no
"the" maximal u-RAF.

```
Food: a, b
r1 : a + b [a] {d} => c      # inhibited by d
r2 : a + b [b] {c} => d      # inhibited by c
```

Two maximal u-RAFs, `{r1}` and `{r2}`: each is an RAF whose support avoids its own
inhibitor, and their union is an RAF that fails the uninhibited condition.

`simulate` respects inhibition too — an inhibited reaction has propensity zero, so a
running network can **lose** a subRAF, not merely gain one. See
`examples/inhibition_dissolution.py`.

The set-theoretic tools need no special handling: `sample_irrraf`, `irrraf_census`,
`core_raf` and `catalytically_reachable` all take a reaction set, and passing a u-RAF
is correct because the uninhibited property is inherited downward — every sub-RAF of a
u-RAF is a u-RAF.

Deciding whether a u-RAF exists is NP-complete, but the problem is fixed-parameter
tractable in *k*, the number of inhibition classes — and **k is a property of how you
encode inhibition, not of the chemistry.** `classes_from_inhibitors` groups by
inhibiting *molecule*, so *k* is the number of distinct inhibitors rather than the
number of inhibited reactions, which is the difference between `2^k` being feasible
and not.

## Notes on irreducible cores

There may be **exponentially many** irreducible RAFs inside one maximal RAF, and
finding the smallest is NP-hard (Steel, Hordijk & Smith 2012). `sample_irrraf` returns
*one*, chosen by the random order it walks; `irrraf_census` samples repeatedly and
reports how many distinct ones it saw. That count is always a **lower bound**, never
an upper one.

`is_food_catalysed` exists because a core whose every reaction has a food catalyst
satisfies the letter of the RAF definition while being in no sense self-referential —
it runs wherever the food runs. Split those out before reading a count of cores as a
count of anything biological.

## Install

```bash
pip install rafkit
```

Development:

```bash
pip install -e ".[test]"
pytest -q
```

Releases are documented in [CHANGELOG.md](CHANGELOG.md); the release procedure is
[docs/RELEASING.md](docs/RELEASING.md).

## Citing

Cite the concept DOI [10.5281/zenodo.21954795](https://doi.org/10.5281/zenodo.21954795),
which always resolves to the latest version; `CITATION.cff` also lists the per-version
DOI. If you use the CatReNet interoperability or the validation fixture, please cite
CatReNet too.

## References

- Hordijk & Steel, "Detecting autocatalytic, self-sustaining sets in chemical reaction
  systems," *J. Theor. Biol.* 227, 451 (2004).
- Steel, Hordijk & Smith, "Minimal autocatalytic networks," *J. Theor. Biol.* 332, 96
  (2013); arXiv:1212.4450.
- Hordijk & Steel, "Autocatalytic sets extended: dynamics, inhibition, and a
  generalization," *J. Syst. Chem.* 3, 5 (2012); arXiv:1206.1017.
- Huson, Xavier & Steel, "CatReNet: interactive analysis of (auto-)catalytic reaction
  networks," *Bioinformatics* 40(8), btae515 (2024).
- Serra & Villani, "Template-Based Catalysis and the Emergence of Collectively
  Autocatalytic Systems," *Entropy* 28(2), 184 (2026).
- Matsubara, Ameta, Thutupalli, Nghe & Krishna, "Conditions for Darwinian evolution in
  compartmentalized autocatalytic reaction networks," arXiv:2211.03155 — the analytic
  benchmark behind `rafkit.dilution`.

## License

MIT. CatReNet is GPL v3 and is **not** a dependency — this library interoperates with
it only through files, and contains no code derived from it.
