Python · chapter 11 of 13

Accuracy and validation

"Accurate generator" is a claim that should be backed by numbers. GRNLab is checked at three levels.

1. Unit tests

tests/ holds 50 pytest tests (pytest in the python folder). They cover analytic steady states of one- and two-stage models, finite-difference checks of every Jacobian, the textbook behaviour of the motifs (bistable toggle, oscillating repressilator and Goodwin loop, pulsing incoherent feed-forward loop, monotone coherent one), seeding, every noise model, perturbations, stimuli, steady states, screens, export round trips, and figure rendering.

2. The accuracy audit

python -m grnlab.audit exercises the public API end to end against closed-form results and writes AUDIT.md. It is independent of the unit tests and meant to be rerun after any change to the kinetics or integrators.

Check	Analytic result	Tolerance
Unregulated gene trajectory	`x(t) = x* + (x0 - x*) e^{-gamma t}`	1e-6
Hill function at threshold	`h(K) = 1/2`, `h'(K) = n / (4K)`	1e-14
Activation-edge steady state	closed form from `x1* = vmax/gamma`	1e-8 relative
Repressilator Hopf threshold (Hill 4)	`alpha_c = 2` (loop slope `n x^n / (1 + x^n) = 2` at `x* = 1`)	1e-6
Below/above threshold	decays / sustained oscillation
Toggle switch bistability threshold (Hill 2)	`alpha_c = 2`	1e-6
Jacobian vs central differences	0	1e-6
Euler–Maruyama convergence	error ratio 2 when `dt` halves	0.3
Additive noise, linear gene	stationary variance `sigma^2 / (2 gamma)`, mean `vmax/gamma`	10 %
Chemical Langevin, linear gene	Fano factor `sigma^2` (linear-noise approximation)	8 %
Two-stage protein level	`q* = beta m* / gammaP`	1e-8
Fast-protein limit	two-stage steady state equals one-stage	1e-5
Step stimulus	exact exponential after the switch, exactly zero before	1e-6
Knockout screen	every knockout with a resolvable effect (above 1e-3 log2 units) moves its target in the sign predicted by `S`	exact
Seeds	two runs are bit-identical	exact

The last run passed 19 of 19 checks; the measured values are in python/AUDIT.md. The MATLAB toolbox runs the same audit and passes it too.

python -m grnlab.audit

Output

check                                                           expected      measured       error       tol  pass
Unregulated gene x(t): max |ODE - analytic|                            0   1.01804e-07    1.02e-07     1e-06  True
Hill h(K) = 1/2                                                      0.5           0.5           0     1e-15  True
Hill dh/dx at K = n/(4K)                                        0.441176      0.441176    5.55e-17     1e-14  True
Activation edge x2* (additive)                                   1.70927       1.70927    2.22e-16  1.71e-08  True
Repressilator Hopf alpha_c (hill 4)                                    2             2    4.82e-14     1e-06  True
Repressilator amplitude at alpha = 1.8 (below threshold)               0   1.97282e-07    1.97e-07     1e-06  True
Repressilator oscillates at alpha = 2.4 (above threshold)              1             1           0         0  True
Toggle switch bistability alpha_c (hill 2)                             2             2    8.88e-16     1e-06  True
Jacobian vs central differences (max abs)                              0   2.58534e-10    2.59e-10     1e-06  True
Euler error ratio for dt halved (expect 2)                             2       2.00065    0.000655       0.3  True
OU stationary variance sigma^2/(2 gamma)                           0.125      0.128342     0.00334     0.012  True
OU stationary mean vmax/gamma                                         20       19.9801      0.0199      0.02  True
Langevin Fano factor var/mean = sigma^2                             0.16      0.161622     0.00162     0.012  True
Two-stage protein q* = beta m*/gammaP                            14.7902       14.7902    1.78e-15     1e-08  True
Fast-protein limit: max |two-stage - one-stage| steady state           0   1.77636e-15    1.78e-15     1e-05  True
Step stimulus x(t): max |ODE - analytic|                               0   2.86285e-08    2.86e-08     1e-06  True
KO effect sign agreement with S (fraction)                             1             1           0         0  True
KO screen convergence (fraction)                                       1             1           0         0  True
Seed reproducibility: max |run A - run B|                              0             0           0         0  True
19 of 19 checks passed

Two of the audited thresholds. Left: repressilator amplitude against alpha; oscillation starts at the analytic Hopf point alpha = 2. Right: toggle-switch steady states against alpha; the symmetric state (blue) loses stability at alpha = 2 and two branches (teal) appear.

Euler–Maruyama with zero noise against RK45: the error falls linearly with the step, as expected for a first-order scheme.

3. Numerical methods and their limits

ODE integration. SciPy's RK45 or BDF with rtol=1e-7, atol=1e-9; solver output is clipped at zero. The integration restarts at every stimulus edge. Sampling is exact at the requested times (the solver is asked for them), not interpolated afterwards.
Steady states. Integration until the relative residual is below tol, then Newton polishing to machine precision. Oscillators are reported as converged=False with a time average.
SDE integration. Euler–Maruyama, strong order 0.5, with clipping at zero. The clipping introduces a bias of order sigma sqrt(dt) near zero; keep dt small (the default is sample_dt / 20) and, for the linear-noise checks above, expression well above the noise amplitude. Milstein is not needed for additive noise and gives no order improvement for the Langevin form with its state-dependent diffusion unless the derivative is included, so it is not offered.
Chemical Langevin. Valid when molecule numbers are large enough that fluctuations are Gaussian; it is not a Gillespie simulation and will not reproduce bursting or extinction statistics at very low copy numbers.
Hill kinetics. The model is a phenomenological transcription model: no explicit binding kinetics, no cooperativity between different regulators beyond the AND/OR combination, no delays other than the mRNA stage.

Known limitations

Additive observation noise is clipped at zero, which biases low values slightly upward; use sol.x_clean when that matters.
The force layout of grn.plot_network depends on its seed; labels can still overlap on dense networks above about 80 genes.
The Euler–Maruyama loop runs in Python; long noisy simulations of large networks are the slowest case.
The multiplicative-logic Jacobian is O(n^3); use solver="ode45" above a few hundred genes.

← Recipes Function reference →