Metadata-Version: 2.4
Name: rafkit
Version: 0.7.1
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, binary polymer chemistries whose catalysis is drawn, structured or templated, and
interoperability with [CatReNet](https://github.com/husonlab/catrenet) and the Petri net
tools.

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's K-BPM: catalysis by coin flip, with optional cleavage |
| `complementary_polymer` | Serra & Villani's C-BPM: catalysis by active-site complementarity, not by coin flip |
| `firing_disk_polymer` | their firing-disk construction — a chemistry *grown* from a seed, closed under its own production |
| `templated_polymer` / `templated_catalysts` | templated ligation: a species catalyses `a + b -> ab` iff it carries the product's complement — nothing is drawn |
| `degree_preserving_null` / `motif_matched_null` | same reactions, rewired catalysis — the nulls a structured chemistry is compared against |
| `matched_f_random` / `matched_f_cbpm` | the K-BPM and C-BPM at a structured chemistry's catalysis level |
| `catalysis_motifs` | what a chemistry and its nulls are compared on |
| `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 |
| `to_andl` / `write_andl` | ANDL export with mass-action rate constants — runs unchanged in Spike |
| `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.

## Catalysis three ways: drawn, structured, and templated

`binary_polymer` is Kauffman's K-BPM: every string up to `max_len` exists, and whether a
species catalyses a reaction is an independent coin flip at probability `p`. Two more
ensembles keep that reaction set and change the one thing that matters — *which* catalyst
catalyses *which* reaction.

**`complementary_polymer`** is Serra & Villani's C-BPM (*Entropy* 28(2), 184, 2026),
reproduced rather than invented. A catalyst carries an active site — a substring of itself,
3–4 residues long by default — and acts on whatever is complementary to that site. A
K-catalyst's targets are independent draws; a C-catalyst's targets all share one template,
so they are structurally correlated, and Serra & Villani measure the signature of that as a
far higher and far more irregular reactions-per-catalyst distribution (~400 against ~20).
**`firing_disk_polymer`** is their other construction. Rather than enumerating every string
and sprinkling catalysis over the result, it grows the chemistry outward from a small seed,
so a species exists only if some reaction actually makes it: an enumerated chemistry is
full of species nothing can reach, and a grown one is closed under its own production by
construction. ⚠ Food is taken to be the firing disk, an assumption the paper leaves open.

**`templated_polymer`** draws nothing. A species templates the ligation `a + b -> ab`
exactly when it contains the (reverse) complement of the product — or, under
`rule="junction"`, of the `2h` residues spanning the junction — given both reactants reach
the per-side overlap `h`. The catalysis graph is a *function* of the sequence set, which is
the content of a template world and the reason this generator has no `rng` argument. The
cleavage carries its ligation's templates, so a template changes a reversible reaction's
rate and never its equilibrium. `templated_catalysts` is the same rule as a pure function
over whatever strands exist, for a simulator that holds explicit strands.

Sized on the complete `max_len` 7 set, the two rules behave very differently, and the
numbers are pinned as known answers:

| | |
|---|---|
| junction rule, `h` ≤ 2 | **saturates** — every eligible reaction templated by 194 or 46 species, f = 978 and 142 against ~5 for a random chemistry at the RAF threshold: a uniform speed-up with no specificity left |
| product rule, `h` 2 | f = 12.9, with the template count *falling* with product length — 46 for a 4-mer, exactly one for a 7-mer, its own complement |
| either rule, `h` 3 | the **maximal RAF is empty**: templated reactions need reactants of length ≥ `h`, and only untemplated reactions make those from a food set shorter than `h` |

The last row is the one to remember. `max_raf` remains the right question about closure
and the wrong one about whether templating acts — the chemistry is perfectly runnable on an
uncatalysed background.

**A structured chemistry needs something to be compared against**, and "random at the same
f" is a poor choice: at a templated f, a random chemistry spreads its edges over nearly
every reaction (reach ~0.92 against the product rule's 0.61). So each null holds something
different fixed, and `catalysis_motifs` reports what they are compared on — f, reach,
catalysts per reaction by product length, self-catalysed reactions, mutual pairs, and the
maximal RAF's size, all counted on the reversible pair:

- `degree_preserving_null` — a double-edge-swap chain keeping every degree;
  `stratified=True` also keeps the *lengths* of each reaction's templates, which the plain
  shuffle breaks along with the sequence tie.
- `motif_matched_null` — random catalysis with exactly the network's edge count,
  self-catalysed reactions and mutual pairs; `match_lengths=True` plants them at the
  network's own lengths, without which it is a second random chemistry with two counts
  pinned.
- `matched_f_random` / `matched_f_cbpm` — the K-BPM and C-BPM at the network's catalysis
  level, in expectation.

Every one of these returns the same `BinaryPolymerNetwork`, so everything in the table
above runs on them unchanged.

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

### Readable is not runnable: `to_andl`

PNML's `ptnet` grammar has no place for a rate constant, so a PNML file documents a
network without being able to run it. `to_andl` / `write_andl` export the same network to
ANDL, the PetriNuts format shared by Snoopy, Spike and Marcie, with mass-action rate
constants — a complete stochastic Petri net that Spike executes directly. An independently
developed simulator re-running this chemistry from its definition is the entire point.

The semantics are mass action and nothing else, and they were **measured** against Spike
1.6.0rc2 rather than read from its documentation: a catalyst is a consume-and-produce
self-loop, so the propensity scales with catalyst count; alternative catalyst sets are
separate transitions whose propensities sum; and `a + a -> aa` is a weight-2 arc counted
as unordered pairs, `n(n−1)/2`, with the alternatives 19–25 standard errors away.
⚠ Catalysis scaling with catalyst count is *not* what `simulate` does — there any catalyst
present buys the full rate — so this export runs a different system from the in-library
simulator, deliberately; the pair convention, on the other hand, is `simulate`'s exactly.
Rate constants are emitted as named constants, so a Spike `.spc` configuration can override
any single rate without regenerating the file.

Refused rather than silently altered, because an executable file that drops a feature does
not *document* a different system, it **runs** one: inhibition (PNML may annotate it, since
a reader sees the annotation; a simulator would not), and a catalyst that is also a
reactant of the same reaction, whose self-loop merges with the consuming arc into a
weight-2 pre-arc with implementation-dependent combinatorics. The latter means a
`templated_polymer` network, in which a template can be its own reactant, is refused until
someone measures what Spike does with it.

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