knockwire

Recover signed interaction networks from knockout experiments, and report which edges the data cannot determine.

The problem

Systems biologists knock down genes and control engineers nudge setpoints, then want the signed interaction matrix back. Modular Response Analysis (MRA) and its regression variants always return one matrix, even when the experiment design cannot tell several networks apart. Perturb half the nodes and a pseudo-inverse or a lasso penalty silently fills in the missing half. The output looks the same either way, and nothing standard tells you which reported edges the data actually pins down and which are artefacts of the regulariser. In the benchmark below, roughly half of the signs plain MRA reports at 30% coverage are wrong.

How it works

Near a steady state, x' = f(x, p) with Jacobian J. Perturbing node k shifts the steady state by the global response column R_k = -J^{-1} e_k d_k. Here d_k is the direct effect on node k, which is usually unknown. The data therefore fix each row of J only up to scale, so knockwire works with MRA's local response matrix r_ij = -J_ij / J_ii (diagonal -1). With every node perturbed this is r = -diag(R^{-1})^{-1} R^{-1}.

A perturbation of node k does not act on any other node i directly. So every perturbed k != i gives one linear equation for row i:

sum_{j != i} r_ij R_jk = R_ik

1. Estimation. Replicates are averaged to fix the design, then kept as separate equations. They are solved with total least squares (smallest right singular vector of [A | b]), restricted to the row space the design can see.

2. Identifiability. The solutions for row i form an affine set theta0 + span(N_i), where N_i is the SVD null space of that row's design. Edge j -> i is identified when row j of N_i is zero: every consistent network has the same value. Otherwise the edge can take any real value, and knockwire marks it unidentifiable. It can return two networks that reproduce the data exactly but give that edge opposite signs (Inference.witnesses).

Optional priors turn the affine set into a polytope:

• a bound on |r_ij|,
• known signs,
• known absent edges.

A hand-written two-phase simplex (Bland's rule, numpy only) then minimises and maximises every undetermined edge over that polytope. If the range lies in one closed half-line, the edge is sign-identified. Inference.completions samples consistent networks from the polytope by hit-and-run.

3. Noise. Pass the relative measurement noise and knockwire runs a parametric bootstrap. It re-estimates from multiplicatively re-noised data and widens every interval by z = 4 bootstrap standard deviations. A sign is claimed only when the interval excludes zero.

4. Time series. Step responses are fitted with one shared discrete-time model y[t+1] = A y[t] + c_e. Then J = logm(A) / dt, using inverse scaling and squaring (Denman-Beavers square roots plus a Mercator series, no SciPy). The steady-state responses (I - A)^{-1} c_e come out exactly, even from trajectories that stop long before settling, and feed straight into step 2.

5. Experiment design. Which edges a new knockout would pin depends on its response pattern, which is unknown before the experiment is run. For each candidate node, knockwire simulates generic networks on the plausible wiring and counts the undetermined edges that would remain. By default the plausible wiring is every edge not already identified as zero; you can also pass a known wiring diagram. The suggestion is the candidate with the lowest count.

A tiny worked case. Perturb only c in a three-node system and measure R = (a: 0.5, b: 0.25, c: 1.0). Row a gives one equation, 0.25 r_ab + 1.0 r_ac = 0.5: a line, so both edges are unidentifiable. MRA's pseudo-inverse reports r_ac = 0.47, r_ab = 0.12 anyway. Add the prior "b represses a" (r_ab <= 0) and the LP proves r_ac >= 0.5, so c -> a becomes sign-identified positive. Add |r| <= 1 and its range tightens to [0.5, 0.75]. This exact case is a unit test and a CLI test.

Synthetic ground truth comes from a built-in generator of random signed sparse networks with linear or Hill-type dynamics. Steady states are found by RK4 integration polished with Newton steps on the analytic Jacobian.

Install and usage

Python 3.10+; the only runtime dependency is numpy.

git clone <this repo> && cd knockwire
python -m pip install -e ".[dev]"
pytest -q

Generate a synthetic experiment: a 5-node Hill network, 3 of 5 nodes knocked down 5%, three replicates, 5% multiplicative noise. Then infer:

$ knockwire simulate --nodes 5 --density 0.4 --coverage 0.6 --noise 0.05 --replicates 3 --seed 7 --out data.csv
wrote data.csv: 5 nodes, perturbed a,d,e
true signs (row = target, column = source):
    a . 0 0 - +
    b 0 . 0 0 0
    c - - . 0 0
    d 0 0 0 . 0
    e + - 0 - .

$ knockwire infer data.csv --noise 0.05
Signed local response matrix (row = target, column = source)
  + / - : sign determined   0 : identified as zero   ? : sign not determined
      a   b   c   d   e
  a   .   ?   ?   ?   ?
  b   ?   .   ?   ?   ?
  c   -   ?   .   0   0
  d   ?   ?   ?   .   ?
  e   ?   ?   ?   ?   .

edges: 3 identified, 0 sign-identified, 17 unidentifiable

edge         estimate                interval  status
b -> a     -3.564e-18             [-inf, inf]  unidentifiable
c -> a       -0.02472             [-inf, inf]  unidentifiable
d -> a        -0.2677             [-inf, inf]  unidentifiable
e -> a         0.1987             [-inf, inf]  unidentifiable
...
a -> c        -0.2039      [-0.2393, -0.1686]  identified
b -> c     -1.093e-17             [-inf, inf]  unidentifiable
d -> c       0.003809     [-0.01032, 0.01794]  identified
e -> c      -0.004661    [-0.01494, 0.005617]  identified
...
a -> e         0.2063             [-inf, inf]  unidentifiable
d -> e        -0.1792             [-inf, inf]  unidentifiable

The minimum-norm estimates for d -> a, e -> a, a -> e and d -> e happen to have the right signs here. With this design they are still guesses, and knockwire says so. A pseudo-inverse reports them with the same confidence as a -> c.

$ knockwire suggest data.csv --noise 0.05
undetermined edges now: 17 (the wiring model predicts 20.0 for this design)
suggested next perturbation: b
  perturb    b: 16.0 undetermined edges predicted
  perturb    c: 16.0 undetermined edges predicted

The CSV format is a header node,<one label per experiment> followed by one row per node. By default the column labels name the perturbed node; --perturbed a,a,b,... overrides that, one entry per column, and repeats mark replicates. Priors: --bound 2, --sign b>a=- (repeatable). --zero-below 1e-6 treats tiny responses as structural zeros.

Library:

import knockwire as kw

net = kw.random_network(12, density=0.2, kind="hill", rng=0)
exp = kw.simulate(net, targets=[0, 3, 5, 8], strength=-0.05, noise=0.1, replicates=3, rng=1)
result = kw.infer(exp.responses, exp.targets, net.names, noise=0.1)
print(result.count(kw.IDENTIFIED), "identified")
status, sign, lower, upper = result.status, result.sign, result.lower, result.upper  # [target, source]
print(kw.suggest_next_perturbation(result).name)

# Exact, noise-free analysis on a linear network, with a magnitude prior.
lin = kw.random_network(8, density=0.3, rng=2)
clean = kw.simulate(lin, targets=[0, 1, 2, 3], strength=0.1)
exact = kw.classify_identifiability(clean.responses, clean.targets, lin.names, bound=5.0)
low, high = exact.witnesses(4, 0)       # two consistent networks, extreme values of edge a -> e
samples = exact.completions(100, rng=3)  # consistent networks sampled from the feasible set

# Time series: step responses sampled every dt, fitted and mapped to J with logm.
trajs = [lin.trajectory(lin.input_for(k, 0.5), dt=0.2, samples=30) for k in range(8)]
fit = kw.infer_timeseries(trajs, list(range(8)), dt=0.2)
print(abs(fit.jacobian - lin.weights).max())

Results

Run python benchmarks/recovery.py (500 networks; 173 s on an 8-core Apple Silicon laptop, arm64, Python 3.12.13, numpy 2). The seed is fixed, so rerunning prints the same table; the raw output is kept in benchmarks/results.md.

Setup: 500 random signed sparse linear networks, 10-30 nodes, mean in-degree 2. Every network is run at every noise level (multiplicative Gaussian, relative sd 0-20%) and every coverage (30/60/100% of nodes perturbed, random subset). Each perturbation has three replicates, and knockwire uses 50 bootstrap draws.

A sign call claims an edge is positive or negative:

• Plain MRA: |r| > 0.05 (exact inverse at full coverage, minimum-norm least squares otherwise).
• Lasso: row-wise coordinate descent with penalty 0.05 * lambda_max; any non-zero coefficient is a call.
• knockwire: only an interval that excludes zero.

A call is false if the edge is absent or has the opposite sign. F1 is computed over true signed edges. The false-call rate is false calls / all calls, with knockwire's raw count of false calls in parentheses.

noise	coverage	F1 knockwire	F1 MRA	F1 lasso	false calls knockwire	false calls MRA	false calls lasso	identified (knockwire)
0%	30%	0.036	0.388	0.300	0.0000 (0)	0.510	0.697	6.9%
0%	60%	0.219	0.629	0.538	0.0000 (0)	0.378	0.507	25.9%
0%	100%	1.000	1.000	0.885	0.0000 (0)	0.000	0.192	100.0%
5%	30%	0.017	0.391	0.297	0.0000 (0)	0.509	0.703	4.3%
5%	60%	0.047	0.629	0.532	0.0000 (0)	0.380	0.518	10.1%
5%	100%	1.000	0.999	0.870	0.0007 (14)	0.000	0.218	100.0%
10%	30%	0.016	0.391	0.298	0.0000 (0)	0.509	0.704	4.4%
10%	60%	0.048	0.628	0.523	0.0000 (0)	0.379	0.531	10.6%
10%	100%	0.999	0.998	0.842	0.0006 (11)	0.002	0.262	100.0%
20%	30%	0.018	0.388	0.293	0.0000 (0)	0.515	0.718	4.7%
20%	60%	0.047	0.624	0.501	0.0000 (0)	0.389	0.573	10.5%
20%	100%	0.938	0.979	0.765	0.0005 (9)	0.037	0.371	100.0%

What the table shows:

• Noise-free data. knockwire never makes a false call. That is exact by construction and holds in every cell.
• Partial coverage. The baselines get their higher F1 by calling signs the data do not determine, and 38-72% of those calls are wrong. knockwire's F1 is low there because it declines to guess. Most of what it does identify at partial coverage are absent edges (identified as zero), and those do not count toward F1.
• Full coverage with noise. knockwire makes 9-14 false calls out of about 20,000, a rate of 0.05-0.07%, so the rate is not zero. At 5-10% noise, MRA with a 0.05 threshold does slightly better (0.0-0.2%); at 20% noise, MRA's rate rises to 3.7% while knockwire's stays at 0.05%. knockwire's remaining false calls come from the bootstrap: its 4-sd intervals are centred on a noisy estimate, so a few true-zero edges land outside them.
• Identified share. At 60% coverage the identified share drops from 26% without noise to about 10% with any noise; see the first limitation below.

Design notes

The central decision was to separate what the design determines from what the estimate says. The null space depends only on which nodes were perturbed and the zero pattern of the responses. The point estimate is just one element of the solution set. Regression methods merge the two, and a regulariser then picks one element of that set for you. Keeping them apart costs something: at partial coverage knockwire refuses to report signs that MRA gets right about half the time. I took that trade because a wrong edge in a signalling model is expensive to discover later, while a "?" is honest and cheap.

Two smaller choices follow the same logic. First, the sign-identification LP uses a hand-written simplex over the null-space coordinates rather than a projected method. Exact vertex solutions double as witnesses, so every "unidentifiable" verdict ships with two concrete networks that prove it. The tests check both witnesses against the data. Second, noise is handled by bootstrap on top of the exact structural analysis rather than by a statistical rank test. The structural answer does not move when the noise level is misstated. The noise level only widens or narrows intervals.

Limitations

• Identification that relies on exact proportionality is lost under noise. Some edges are pinned only because two nodes respond in exact proportion, for example both sitting downstream of a single bottleneck. Noise breaks that proportion, so knockwire then reports those edges as unidentifiable. The benchmark's "identified" column falls from 26% to 10% at 60% coverage as soon as any noise is added, for this reason. The error is conservative (fewer claims, not wrong ones), but it does lose information.
• Structural zeros must be exact. Zeros in the responses (nodes a perturbation cannot reach) drive most identification at partial coverage. Multiplicative noise preserves them; real assays report small non-zero values. Use --zero-below with a detection threshold, and accept that a bad threshold changes the answer.
• The noise model is assumed, not estimated. You supply the relative noise level. The bootstrap resamples around the noisy estimate, so intervals are approximate: the table shows a small number of false calls at full coverage with noise.
• Linearisation. Finite perturbations of nonlinear systems bias the estimates at second order. With noise=0 on Hill data, absent edges come out around 1e-3 and get a sign; always pass a noise level for real data.
• Sign-identification needs priors. With no bound or known signs, every undetermined edge ranges over the whole real line, so the sign-identified class is empty. The benchmark runs without priors.
• Scale. The LP runs two simplex solves per undetermined edge per row, in pure numpy. That is fine for tens of nodes and slow for hundreds, especially combined with the bootstrap.
• Experiment design is only as good as the assumed wiring. Given the true wiring diagram, the suggestion is within one edge of brute force (tested against brute force on 40 random networks of 3-6 nodes). Without it, knockwire assumes every undetermined edge may exist, and the prediction can be off.
• Time series. The shared linear model assumes step inputs from the unperturbed steady state. The Jacobian is returned only when the trajectories span the state space, and logm requires dt small enough that no eigenvalue of J dt has imaginary part beyond pi.
• The benchmark uses linear networks for speed; Hill networks are covered by the test suite, not the benchmark.

License

MIT. See LICENSE.