MATLAB · 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 60 matlab.unittest tests (buildtool runs them, and buildtool package refuses to package if any fail). 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

audit/runAudit.m exercises the public API end to end against closed-form results and writes audit/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 audit/AUDIT.md. To rerun:

cd audit; report = runAudit();

Output

                                <strong>check</strong>                                 <strong>expected</strong>     <strong>measured</strong>       <strong>error</strong>       <strong>tolerance</strong>    <strong>pass</strong> 
    <strong>______________________________________________________________</strong>    <strong>________</strong>    <strong>__________</strong>    <strong>__________</strong>    <strong>_________</strong>    <strong>_____</strong>

    "Unregulated gene x(t): max |ODE - analytic|"                           0     1.1862e-07    1.1862e-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.44118        0.44118    5.5511e-17        1e-14    true 
    "Activation edge x2* (additive)"                                    2.297          2.297             0    2.297e-08    true 
    "Repressilator Hopf alpha_c (hill 4)"                                   2              2    4.4409e-16        1e-06    true 
    "Repressilator amplitude at alpha = 1.8 (below threshold)"              0     2.2272e-07    2.2272e-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             0        1e-06    true 
    "Jacobian vs central differences (max abs)"                             0     2.7891e-10    2.7891e-10        1e-06    true 
    "Euler error ratio for dt halved (expect 2)"                            2         2.0012     0.0012015          0.3    true 
    "OU stationary variance sigma^2/(2 gamma)"                          0.125        0.13381     0.0088071        0.012    true 
    "OU stationary mean vmax/gamma"                                        20          20.01      0.009728         0.02    true 
    "Langevin Fano factor var/mean = sigma^2"                            0.16        0.15901    0.00099419        0.012    true 
    "Two-stage protein q* = beta m*/gammaP"                            6.4843         6.4843    8.8818e-16        1e-08    true 
    "Fast-protein limit: max |two-stage - one-stage| steady state"          0     4.4409e-16    4.4409e-16        1e-05    true 
    "Step stimulus x(t): max |ODE - analytic|"                              0     1.7995e-08    1.7995e-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

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 ode45: the error falls linearly with the step, as expected for a first-order scheme.

3. Numerical methods and their limits

ODE integration. ode45 or ode15s with RelTol 1e-7, AbsTol 1e-9, and NonNegative on every state. 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 SampleDt / 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.xClean when that matters.
The force layout of grn.plotNetwork depends on its Seed; labels can still overlap on dense networks above about 80 genes.
The multiplicative-logic Jacobian is O(n^3); use ode45 above a few hundred genes.

← Recipes Function reference →