biosym

Evolutionary discovery of interpretable biological models.

biosym searches for compact closed-form equations that explain experimental data — Michaelis–Menten kinetics, Hill dose-response curves, growth laws — using genetic-programming symbolic regression from coevo. It returns a whole Pareto front of models trading off accuracy against complexity, so you can pick the simplest law that still fits.

The idea

In biology the goal is usually understanding, not black-box prediction. A compact equation like V·x / (K + x) tells you the mechanism; a neural net doesn't. biosym searches the space of closed-form expressions, biased by a library of biologically-motivated operators, and shows you every model on the accuracy-vs-simplicity frontier.

Fitting a curve and recovering a law are different achievements, and this README tries to keep them apart. Every accuracy number below is measured on data the model was not fitted to, and there is a separate column for whether the discovered expression is actually the ground-truth law.

Install

pip install -e ".[dev]"          # numpy, coevo, pytest, sympy

discover() needs only numpy and coevo. sympy is used to score symbolic recovery (pip install "biosym[symbolic]"); PySR is optional and only used as a benchmark baseline.

Quickstart

from biosym import datasets, discover

X, y = datasets.hill(n=150, V=1.0, K=2.0, n_hill=2.0, noise=0.02)
front = discover(X, y)

for model in front:
    print(f"rmse={model.rmse:.4f}  complexity={model.complexity}  y = {model.expression}")

# every model on the front can be applied to new data
best = front[-1]
predictions = best.predict(X_new)
held_out_rmse = best.score(X_test, y_test)

python benchmarks/run_discovery.py           # recover synthetic laws
python benchmarks/compare_real.py --pysr     # vs baselines and PySR, held out

What's inside

Module	Purpose
operators	biological operator library (pow, sigmoid, log1p, expm1, relu)
datasets	synthetic ground-truth laws: Hill, Michaelis–Menten, logistic, exponential, power law, damped oscillation
discover	multi-objective discovery → (RMSE, complexity) Pareto front of applicable models
symbolic	is a discovered expression the ground-truth law, or just a good fit?
benchmark	held-out comparison against linear/polynomial baselines and PySR
ode	discovering governing dynamics from time-course data

How discovery works

discover runs coevo.SymbolicRegressor, a genetic-programming engine that evolves expression trees with subtree crossover and mutation. Two mechanics do most of the work:

• Linear scaling (Keijzer, EuroGP 2003) scores a candidate on how well a + b·f(x) fits the target, with a and b obtained in closed form. The search only has to find the shape of the relationship and gets scale and offset for free.
• A per-complexity hall of fame records the best expression found at every tree size, so a single run yields the whole accuracy/complexity front instead of one winner.

Choosing an operator set

The biological operators were intended as a search prior — the idea being that a library of growth/saturation/dose-response shapes would help find biological laws. Measured against the alternative, that turns out to be backwards for the thing this library is for.

Replicated on seeds 20–39 (disjoint from the seeds that generated the hypothesis) across all seven benchmark datasets, 140 runs per arm:

operator set	test RMSE	solved (exact)	solved (form)
plain (default)	0.0916	20/140	45/140 (32%)
biological	0.0800	20/140	23/140 (16%)

The biological operators fit better (Wilcoxon p = 8.8e-06) and recover the ground-truth form half as often (Fisher p = 0.0032). Both clear a Bonferroni-corrected threshold. Exact recovery is unchanged.

The effect concentrates in the exponential-family laws, which makes the mechanism legible:

dataset	biological → plain
exponential_association	1/20 → 12/20
pharmacokinetic_decay	1/20 → 7/20
michaelis_menten	0/20 → 3/20
hill_dose_response	1/20 → 3/20
kleiber_allometry	20/20 → 20/20 (a line in log-log space)
logistic_growth	0/20 → 0/20 (neither recovers it)
gompertz_growth	0/20 → 0/20 (neither recovers it)

expm1 and log1p supply near-miss approximations of an exponential shape, so the search reaches a good fit and never needs the clean exp(-k·t) form. Strip them and the plain exponential is the only route there. The same argument applies to sigmoid and rational saturation.

front = discover(X, y)                            # plain — optimise for mechanism
front = discover(X, y, operators="biological")    # optimise for fit

So "plain" is the default, which is an awkward thing for a library called biosym to conclude — but this README's opening argues that mechanism is the point, and the biological operators trade mechanism for fit. The operator library is still there, still supported, and still the right choice when a close fit matters more than the form.

Recovering known biological laws

Against noisy synthetic data (noise=0.02), five of the six laws are recovered to within the noise floor. The sixth is not, and is reported here rather than omitted — the failures are the informative rows.

ground truth	RMSE	at noise level?
Hill V·xⁿ/(Kⁿ+xⁿ) (n=2)	0.0226	yes
Michaelis–Menten V·x/(K+x)	0.0212	yes
logistic growth L/(1+e⁻ᵏᵗ)	0.0212	yes
exponential decay A·e⁻ᵏᵗ	0.0213	yes
allometric power law a·xᵇ	0.0234	yes
damped oscillation A·e⁻ᵞᵗ·sin(ωt)	0.1630	no — 8× noise

Damped oscillation needs a frequency constant (ω=3) inside a periodic function, which is the classic hard case for genetic programming: the fitness landscape in ω is oscillatory, so no local search can walk to it and the constant has to be found by luck.

run_discovery.py prints the recovered expressions. They are worth looking at, because reaching the noise floor does not mean the expression is the law — the power-law row fits to within noise using 41 nodes of machinery that happens to trace the right curve.

Benchmark: held out, over seeds, against PySR

benchmarks/compare_real.py fits every method on 75% of each dataset and scores it on the remaining 25%, over five data seeds, on five literature-derived biological datasets. Targets are z-scored so RMSE is comparable across datasets.

method	test RMSE	train RMSE	overfit	solved (exact)	solved (form)	secs
pysr	0.0670	0.0576	1.16×	20%	40%	3.5
biosym	0.0698	0.0584	1.19×	12%	12%	4.1
poly-3	0.1141	0.1052	1.09×	—	—	0.0
linear	0.4991	0.4719	1.06×	—	—	0.0

The honest picture. On held-out accuracy biosym and PySR are a dead heat — 12 wins, 13 losses, 0 ties over 25 paired runs, Wilcoxon p = 0.426 — with biosym marginally behind on the median and about 1.2× slower. Both comfortably beat a cubic polynomial.

On recovering the law, PySR is clearly ahead, and that is the metric this library exists to move. It finds the right functional form more than three times as often.

The most useful detail is that biosym's two solution rates are equal while PySR's doubles. PySR's near-misses are the right mechanism with imprecise parameters — it recovers Michaelis–Menten as 1.28 − 1.765/(x + 0.555) against a true K of 0.6. biosym's near-misses are structurally wrong, so relaxing the constants buys it nothing. That points at optimising constants during evolution rather than only after it, which is the top of the roadmap.

Both rates are low in absolute terms. On these five datasets, neither method usually recovers the biological law — it fits it.

Reproducibility note: PySR runs with deterministic=True, parallelism="serial", because random_state alone does not make its search reproducible. That costs it a little accuracy (0.0631 → 0.0670) and about 4× wall-clock. PySR's own default model_selection="best" is used; its "score" setting is far more parsimony-hungry than biosym's knee rule and would handicap it rather than match it.

Discovering dynamics (evolutionary SINDy)

biosym.ode.discover_ode estimates derivatives from a time series and symbolically regresses dx/dt against the state:

from biosym import datasets
from biosym.ode import discover_ode

t, X = datasets.logistic_ode()        # dx/dt = r*x*(1 - x/K)
system = discover_ode(X, t)
print(system["x0"].expression)        # -> (x0 - (x0)**2)

system	discovered dynamics
exponential decay dx/dt=-kx	dx0/dt = -(x0·sigmoid(0))
logistic dx/dt = r·x·(1-x/K)	dx0/dt = (x0 − x0²)
harmonic oscillator	dx0/dt = x1, dx1/dt = −x0

This works only on clean trajectories. Derivatives are estimated with central differences, which amplify noise. Measured on the logistic system:

measurement noise (relative)	discovered dx₀/dt
0.000	x0 − x0² ✓
0.003	x0 − x0² ✓
0.017	x0 − x0² ✓
0.067	0.09145 — collapsed to a constant
0.168	0.09402 — collapsed to a constant

Real biological time courses are noisier than the point where this breaks, so treat the module as a demonstration rather than a tool until the weak-form estimator on the roadmap lands.

Roadmap

Ordered by expected effect on the solution rate, which is the number that matters.

• Find out why logistic and Gompertz growth are never recovered — 0/20 under both operator sets, while the exponential family responds strongly to the operator choice. Two laws that no configuration touches is a sharper lead than a percentage point of RMSE.
• Optimise constants during evolution — built and measured in coevo; it moved the solution rate not at all and ships disabled. Linear scaling already supplies the constants that mattered.
• Weak-form / integral derivative estimation for ode, which avoids differentiating noisy data at all (WSINDy). Bigger jump than smoothing.
• Dimensional analysis and monotonicity as hard constraints, pruning the large fraction of the search space that is biologically meaningless.
• Benchmarks vs. SISSO on biological tasks.

Already done (in coevo): linear scaling, per-complexity hall of fame, algebraic simplification, multi-restart constant refinement, size-bounded crossover, and pluggable/removable operator sets.

License

MIT — see LICENSE.