nugap 0.3.0

Vinnicombe nu-gap metric and a pipeline for comparing time-course data across two conditions

nugap

A Python implementation of the Vinnicombe nu-gap metric (δν) and a pipeline that uses it to find which variables change their dynamics between two experimental conditions, from time-course data.

There is no nu-gap implementation in the standard Python control ecosystem — it has lived almost exclusively in MATLAB's Robust Control Toolbox (gapmetric). This package provides a tested, dependency-light one (numpy + scipy + pandas only).

What the nu-gap metric is

For two linear systems P1 and P2, δν(P1, P2) is a number in [0, 1]:

• 0 — identical dynamics,
• near 1 — very different dynamics.

It is computed from the chordal distance between the two frequency responses, gated by a winding-number (topological) condition. Unlike a naive comparison of fitted parameters, it is a true metric, it is bounded, and it correctly treats systems that look very different on paper but behave similarly (and vice versa).

Reference: G. Vinnicombe, Frequency domain uncertainty and the graph topology, IEEE TAC 38 (1993) 1371–1383.

The idea for your application

You have time-course data for thousands of variables under two conditions. The analysis is a pairwise dynamic network: within each condition, every variable is treated as a candidate input for every other variable, and a first-order input->output model is fitted for each ordered pair (i -> j). Then the nu-gap compares condition A's model with condition B's model for each edge. Edges with a large nu-gap are interactions whose dynamics changed.

from nugap import compare_network

# data_A, data_B: dict  variable name -> array (n_replicates x n_timepoints)
edges = compare_network(
    data_A, data_B, t,
    order=1,          # model poles: 1 = first-order, 2 = second-order
    n_zeros=None,     # numerator zeros (default 0 -> all-pole model)
    n=256,            # contour resolution; 256 is plenty for low order
    min_r2=0.5,       # only test pairs with a real relationship
    gate="either",    # how to combine the two conditions' fit quality
    global_null=True, # pool within-condition nu-gaps -> p_global, q_global
)
# one row per edge (source, target, nu_gap, within_median, separation,
# max_r2, q_global), sorted by significance. Flag changes with q_global < 0.1.

order / n_zeros choose the per-edge model. order=1, n_zeros=0 (default) is the first-order K/(τs+1); order=2 is a two-pole system, with n_zeros=1 if you want a zero (the discrete analogue of MATLAB tfest(data, 2, 1)). Higher order needs more time points per trajectory.

gate controls the fit-quality gate across the two conditions: "either" keeps an edge if the relationship is well fit in at least one condition (so relationships that appear or disappear are tested); "both" requires a good fit in both conditions; "mean" uses the mean R² over all replicate fits.

Two things that matter a lot here

1. Gate on fit quality (min_r2). Most variable pairs have no real first-order relationship; those fits are meaningless and produce large, high-variance nu-gaps that would swamp the null. compare_network only tests an edge if a first-order relationship actually holds (R^2 above min_r2) in at least one condition. This is essential — without it nothing is detectable.

2. Confounding. Pairwise first-order identification assumes the i->j relationship is approximately self-contained. In a densely coupled system each output depends on many inputs, so a single pairwise model is misspecified and the within-condition noise floor rises. This is a property of the method (the same one you ran in MATLAB), not the metric. Sparse or modular systems behave well; dense ones need care.

Scale

N variables -> N*(N-1) ordered edges (a million at N=1000). Each fit is a 2-parameter least squares; the cost is the metric, kept cheap by the small n. The per-edge work is independent, so wrap the edge loop in joblib.Parallel / multiprocessing for real datasets, and/or pass include_pairs= to test only a prescreened candidate set.

Install

pip install -e .

Usage

The metric on two systems directly:

from nugap import tf, nu_gap

P1 = tf([1.0], [1.0, 1.0])      # 1/(s+1), continuous
P2 = tf([1.0], [1.0, 1.2])      # 1/(s+1.2)
print(nu_gap(P1, P2))           # ~0.07

# discrete systems use dt; the metric uses the unit circle automatically
Pd = tf([0.5], [1.0, -0.5], dt=0.1)

The full two-condition comparison:

from nugap import compare_conditions

# data_A, data_B: dict mapping variable name -> trajectory (sampled at t)
# u_A, u_B: the known stimulus, if you have one (else omit -> Prony fit)
df = compare_conditions(
    data_A, data_B, t,
    u_A=u, u_B=u,           # drop these for output-only data
    orders=range(1, 5),     # candidate model orders (AIC-selected)
    method="arx",           # or "prony" (output-only), or "auto"
    min_r2=0.9,             # flag variables with poor fits
)
# df is sorted by nu_gap descending, with fit quality per condition

See examples/demo_network.py for a small validated example, and examples/demo_clock.py for a full biological test case: synthetic circadian RNA-seq for 50 genes (core clock genes, clock outputs, CLOCK-independent rhythmic genes, and background) under wild-type vs CLOCK knockout. The KO collapses the cell-autonomous oscillation (a Hopf bifurcation to a damped fixed point) while CLOCK mRNA is still expressed; compare_network recovers the collapse, flagging ~90% of clock/output edges as changed while leaving the CLOCK-independent rhythmic edges (rhythmic in both conditions) and background alone. The generator is examples/clock_sim.py.

The null and short/flat conditions

null_from_reliable_only=True (default) builds the within-condition null only from edges where a real relationship exists in that condition. This matters whenever one condition loses dynamics (e.g. genes go flat in a knockout): fitting noise-to-noise there produces large, meaningless within-condition nu-gaps that would otherwise inflate the null and hide the real changes.

Visualising results

nugap.viz (needs matplotlib + networkx: pip install nugap[viz]) provides three views of a compare_network edge table:

from nugap.viz import volcano, hub_barplot, hub_network

volcano(df, q_thresh=0.1)        # effect size (nu_gap) vs -log10 FDR; best overview
hub_barplot(df, top=20)          # genes ranked by # of significant changed edges
hub_network(df, top_hubs=15,     # directed graph of the most-rewired genes;
            node_groups=classes) # nodes sized by degree, edges coloured by nu_gap

The volcano is the recommended default (scale-independent, shows everything at once). hub_barplot gives the node-level summary that is usually the most interpretable. hub_network draws the most-rewired genes and the significant changed edges among them as a directed graph (pass node_groups, a dict gene -> label, to colour by class); top_hubs controls how many genes appear, so you can keep it readable on dense results. changed_edge_counts(df) returns the per-gene counts behind both hub views. Each plotting function returns a matplotlib Axes so you can compose or restyle. examples/demo_clock.py writes clock_viz.png (volcano + hubs) and clock_hub_network.png.

Modules

• nugap.metric — the nu-gap metric, chordal distance, winding condition.
• nugap.systems — lightweight SISO LTI type (tf, from_zpk, from_control).
• nugap.fitting — identify discrete LTI models from data (ARX / Prony, with AIC order selection). Swap this out to match your MATLAB procedure.
• nugap.network — pairwise (input->output) network comparison across conditions with fit-quality gating and FDR; the main entry point for your application. compare_network, fit_first_order.
• nugap.pipeline / nugap.replicates — single-variable comparison (one model per variable), with and without replicates.

Comparing conditions with replicates (recommended)

If you have replicates, use the replicate-aware pipeline. It fits a model to every replicate, then uses within-condition nu-gaps (replicate vs replicate of the same condition) as a noise floor and compares the between-condition nu-gap against it:

from nugap import compare_conditions_replicates

# reps_A[var], reps_B[var]: 2D array (n_replicates x n_timepoints)
df = compare_conditions_replicates(
    reps_A, reps_B, t, u=u,
    orders=[1],          # fix the order low for short series (see below)
    method="arx",
    global_null=True,    # pool within-condition gaps across all variables
)
# columns include between_median, within_median, p_global, q_global (BH-FDR)
# sorted by q_global ascending; flag changes with e.g. q_global < 0.1

global_null=True pools the within-condition noise across all of your variables into one well-estimated null, which is far more powerful than the ~3 within-pairs a single variable provides.

The most important knob: model order

With few time points, fix the model order low (1, sometimes 2) rather than letting AIC roam. A too-high order makes the per-replicate fit unstable, so replicates of identical dynamics produce a large nu-gap — that variance lands directly in your noise floor and destroys sensitivity. In the bundled demo, order 2 on 14 points gives a within-condition median nu-gap of ~0.16 and detects nothing; order 1 gives ~0.04 and recovers every true change at FDR < 0.1.

Diagnostic: look at within_median / the within-condition null. If it is large (say > ~0.1), your fits are too unstable — lower the order, average replicates, or get more points before trusting the between-condition results.

Matching your MATLAB tfest workflow

tfest needs an input and a response, and returns a continuous transfer function; you then ran gapmetric on those. Here, fit with method="arx" (input/output) and, if you want continuous-domain numbers to match MATLAB, map each fitted model with nugap.to_continuous before comparing. The ranking of variables is essentially the same in discrete or continuous form, so for discovery you can stay discrete.

Known limitations (v0.1)

• SISO only (one signal per variable). MIMO would need the determinant form of the winding condition.
• Systems with poles exactly on the stability boundary (pure integrators / undamped oscillators) are an edge case in the winding condition; fitted models from real data essentially never hit this.
• The fitting layer is intentionally basic. For best results, match the model class and order you used in MATLAB.

License

MIT — see the LICENSE file. (Fill in the copyright holder in that file before publishing.)

Frequency-band and DC-gain options (v0.2)

nu_gap can restrict the metric to an oscillation frequency band and skip the winding-number test, reproducing the frequency-local nu-gap of DyDE (Mombaerts et al. 2019):

from numpy import pi
nu_gap(P1, P2, band=(2*pi/36, 2*pi/16), check_winding=False)  # periods 16-36 h

The default (band=None, check_winding=True) is the full-spectrum Vinnicombe nu-gap with the winding-number condition. fit_first_order(..., min_dc_gain=0.1) rejects models whose absolute DC gain is at/below the floor (DyDE's low-DC-gain filter); dc_gain(model) returns the static input->output gain.

Numerical validation

Correctness of the ν-gap is checked by four mutually independent routes (validation/validate_numerics.py), so the maths is tested rather than a formula re-checked against itself:

1. Closed-form ground truth — exact values computable by hand (constant-gain pairs; the stable-vs-unstable winding case that must give exactly 1).
2. Algebraic invariance — the chordal metric is invariant under w → 1/w, so nu_gap(P1,P2) == nu_gap(1/P1,1/P2) for stable minimum-phase systems (magnitude-independent; shares no code with the package).
3. Independent reference implementation — a from-scratch evaluation (own contour, winding number and pole counts) agrees with the package to < 1e-6 over hundreds of random systems.
4. Robust-stability theorem — the defining property (Vinnicombe): for every controller stabilising both plants, arcsin b(P2,C) ≥ arcsin b(P1,C) − arcsin δν(P1,P2), with the stability margin b(P,C) computed independently from the gang-of-four H∞ norm. This validates the meaning of the number.

The harness includes a --mutate mutation test that perturbs the metric to confirm the checks actually fail on a wrong answer. A fast subset runs in CI (tests/test_numerics.py).