SPIDER -- Stitched Power-spectra for Inferring Directed information flow from incomplete and asynchronous Experimental Recordings

TL;DR: No experiment records the whole brain at once — Neuropixels sessions, calcium imaging fields of view, and iEEG implants each cover a different sliver, with no shared clock across sessions, animals, or labs. SPIDER shows that you can nonetheless recover directed, frequency-resolved information flow across the whole network, because the cross-spectrum of any pair of channels is a marginal, alignment-free quantity: estimate it wherever the pair happens to be co-recorded, stitch all pairs into one global spectral matrix, fill never-co-observed entries with nuclear-norm completion, and then extract direction via minimum-phase spectral factorization and partial directed coherence. Validated on simulations, calcium imaging, 43 IBL Neuropixels sessions from 12 labs (50 areas, anatomical AUC ≈ 0.9), and 43 human iEEG patients, the method surfaces a striking result: spontaneous brain-wide flow is mostly recurrent, except in the theta band, where it organizes into a feedforward hierarchy with the hippocampal formation at its source — in both mouse and human.

The problem: the brain is only ever seen through a keyhole

Directed ("effective") connectivity methods — Granger causality, dynamic causal modeling, partial directed coherence — all assume you recorded every relevant region simultaneously on a common clock. That assumption is dead on arrival for modern large-scale data. A Neuropixels insertion hits a handful of areas; the next session, in a different mouse in a different lab, hits a different handful. Epilepsy patients get electrodes wherever their clinical picture dictates. During spontaneous activity there is no stimulus to align to, so recordings from different sessions cannot even be put on the same time axis.

Worse, partial observation isn't just an efficiency loss — it actively corrupts causal inference. If a hidden region X_3 drives both X_1 and X_2, a recording that sees only X_1 and X_2 will confidently infer a direct X_1 \leftrightarrow X_2 interaction that does not exist. This latent-confounder problem is why you can't just run PDC session-by-session and average the results: each session estimates a different, wrong quantity.

Existing workarounds fall into two camps, both inadequate. Parametric latent-variable models (VAR/GLM stitching) impose rigid dynamics and don't give frequency-resolved answers — a real loss in neuroscience, where theta and gamma interactions carry different functional signals. Covariance "graph quilting" methods can merge partial covariance matrices, but covariances are symmetric: direction is gone.

The key insight: direction lives in the spectrum, and spectra don't need a clock

Here's the observation everything hinges on. For a wide-sense stationary process, two facts hold simultaneously:

1. The cross-spectrum S_{ij}(\omega) between channels i and j is a marginal property of that pair. It doesn't depend on which other channels were in the recording, and — because it's a second-order statistic of a stationary process — it doesn't depend on when the recording happened. Any block that contains both i and j gives an unbiased estimate of the same quantity. No alignment needed across blocks; you only need a valid time axis within each block.

2. Direction is a global property of the assembled matrix, recoverable from it alone. For a stationary, purely non-deterministic process whose spectral matrix S(\omega) is bounded away from zero, the Wiener–Hopf theory guarantees a unique minimum-phase factorization S(\omega) = H(\omega)\,\Sigma\,H(\omega)^H, where H(\omega) is the transfer function and \Sigma the innovation covariance. Inverting gives \bar{A}(\omega) = H(\omega)^{-1}, exactly the frequency-domain VAR transfer function — without ever fitting a VAR. From this you compute partial directed coherence:

the direct, frequency-specific influence of series \ell on series k after conditioning on everything else. For Gaussian processes, integrating -\tfrac{1}{4\pi}\int \log(1-|\pi_{k\leftarrow\ell}|^2)\,d\omega gives a directed information flow in nats.

So the recipe is: aggregate the marginal quantities (pairwise cross-spectra) across blocks, then evaluate the global functional (factorization + PDC) once on the assembled matrix. The asynchrony that kills time-domain methods simply doesn't enter, because phase relationships between i and j — which is where directional information hides — are measured locally within each block.

The pipeline, concretely

Step 1 — Stitch. For each pair (i,j), average the multitaper cross-spectral estimates over all blocks containing both channels. For spike trains, the ordinary PSD is replaced by the Bartlett spectrum (second-order statistics of the firing intensity), so the same machinery covers continuous signals, point processes (including Hawkes-type spiking models), and mixtures.

Step 2 — Complete (if needed). When some pairs were never co-recorded, the matrix has holes. Completion is only well-posed with a structural prior; SPIDER uses a near-low-rank one — population activity driven by a few shared latent sources, S(\omega) = A S_Z(\omega) A^H + S_\varepsilon(\omega) — and solves a nuclear-norm minimization frequency-by-frequency on a real-symmetric augmentation of the Hermitian matrix. Proposition 3 gives checkable sufficient conditions: spectral energy spread across pairs ("non-spikiness") and coverage |\mathcal{O}| \gtrsim rK\log K \cdot \log T, where r is the effective rank and K the number of areas.

Step 3 — Regularize and factorize. Finite-sample stitched matrices need not be positive definite, so a complex Hermitian graphical lasso (\lambda = 0.01 in practice) regularizes the inverse spectrum; Wilson's algorithm then computes the canonical factorization, and PDC drops out. Everything is O(K^3) per frequency and embarrassingly parallel across frequencies.

Two engineering moves make this work on real electrophysiology. First, since the neurons differ across sessions but the areas recur, each area is compressed to one canonical population signal via spectral PCA (the leading eigenvector of the within-area cross-spectral matrix at each frequency), and stitching happens at the area level with coherence-normalized cross-spectra. Second, sessions whose spectral profiles disagree with a reference are dropped — the empirical enforcement of the theory's key assumption, that all blocks sample approximately the same stationary process. Theorems 1 and 2 give uniform consistency, notably in regimes where K grows nearly exponentially in the per-block sample length.

How well does it work?

On a three-node VAR system with a hidden-confounder motif, SPIDER matches the fully-simultaneous estimate almost exactly, while the naive per-block baseline invents the spurious edge. On a nonlinear Izhikevich spiking network — outside the linear theory — SPIDER's recovered flow correlates with ground-truth synaptic influence at r = 0.96, actually beating the full-simultaneous nonparametric estimate (r = 0.89). In 50-dimensional simulations, completion beats zero-filling up to roughly 70% missingness, beyond which the prior can't compensate and performance collapses — a boundary that matches the coverage condition and that the authors then respect in the real-data analyses.

The real-data validations tell one consistent story: observed pairs good, imputed pairs marginal.

The headline validation is the IBL Neuropixels analysis: 43 sessions, 12 laboratories, 39 mice, 4,074 neurons, 50 brain areas — never recorded together — stitched into one directed network. Compared against the Allen anterograde tracing atlas, the strongest inferred flows match anatomical projections with ROC-AUC approaching 0.9. Crucially, this correspondence is carried entirely by directly co-observed pairs; imputed pairs sit near chance. The authors say this plainly, which is to their credit — and it defines exactly how much to trust the completion step (answer: it stabilizes the global estimate and partially recovers structure, but individual imputed edges are not reliable findings).

The scientific payoff: a theta-band feedforward hierarchy

With a validated brain-wide directed network in hand, the authors ask a question no single session could answer: is spontaneous whole-brain flow hierarchical or recurrent? They quantify this with trophic incoherence F_0 \in [0,1] — assign each area a "level" from the flow, then measure how far edges deviate from a perfect one-step feedforward ladder (F_0 = 0: perfect hierarchy; F_0 \to 1: directionless), tested against a null that scrambles edge directions while preserving undirected strengths.

Integrated over all frequencies, the mouse network is essentially non-hierarchical (F_0 \approx 0.92, p \approx 0.10). But resolved by frequency, a single sharp minimum appears in the theta band (4–8 Hz): F_0 = 0.72, z = -2.4, p = 0.001, Bonferroni-surviving, and robust to leave-one-area-out jackknifing (Spearman \rho = 0.99). The polarity is anatomically coherent: subiculum, dentate gyrus, and CA1 — the hippocampal formation — sit at the source end, with midline thalamus, periaqueductal gray, and infralimbic cortex as sinks. Decomposing the flow, the feedforward share peaks in theta (about half the total, with 92% of theta net flow ascending the hierarchy) and collapses toward gamma, where reciprocal flow reaches 84%.

Then the cross-species replication: the same pipeline on resting human iEEG (43 patients, 38 left-hemisphere AAL regions, 71% of pairs directly co-observed across patients who were never recorded together) yields the same signature — theta strongly hierarchical (F_0 = 0.47, z = -15.5), gamma nearly recurrent (F_0 = 0.87). The hippocampus itself is excluded (subcortical), but the cortical sources are exactly its paralimbic partners: olfactory cortex, temporal pole, insula, medial/orbital prefrontal cortex.

This is a genuinely emergent result: frequency-specific directionality — feedforward in theta, recurrent in gamma — is invisible to time-domain or scalar connectivity measures, and invisible to any single recording. It fits the established picture of theta coordinating hippocampal–cortical communication, but here it is read out from the data with no anatomical prior, across species and modalities (spikes vs. field potentials).

What to be skeptical about

The estimand is linear-Gaussian directed structure. PDC via minimum-phase factorization captures everything a second-order-stationary description can. Nonlinear or cross-frequency directed interactions that leave the cross-spectral phase untouched are invisible. The Izhikevich result suggests this is often good enough in practice, but it's a three-node toy.

Exchangeability across animals and patients is a strong assumption. Stitching treats every session as a sample from the same stationary process. Brain states, behavioral context, anesthesia depth, and disease (the iEEG patients have epilepsy, and electrode placement is pathology-driven) all violate this to some degree. The session-selection step filters gross violations but also risks circularity — sessions are kept partly because they agree.

Imputed edges are weak, and the authors' framing depends on that not mattering. The theta-hierarchy analysis is restricted to co-observed pairs, which is the right call. But the human high-band results, where a weak residual hierarchy cleared the null, are flagged by the authors themselves as possibly a completion artifact — the triage question of whether matrix completion can inject spurious directed flow is live, and the near-chance anatomical AUC for imputed pairs says yes, at the edge level it can. Treat any conclusion resting on never-co-observed pairs as provisional.

The mouse theta effect is real but not enormous (F_0 = 0.72 vs. null, z = -2.4) — the human replication (z = -15.5) is what makes it convincing. Also note the PAG-as-apex-sink observation rests on only six co-observed partners; the authors appropriately soft-pedal it.

No comparison to parametric latent-variable stitching methods on shared data; the baselines are zero-fill and per-block PDC, both straw-ish. A head-to-head with, e.g., latent VAR Granger approaches would sharpen the case.

Why this matters

If it holds up, SPIDER converts a structural fact of neuroscience — that brain-wide initiatives produce fragmented, asynchronous recordings by necessity — from an obstacle into raw material. The IBL dataset already existed; nobody could do directed-flow inference on it across sessions until now. The recipe generalizes to any wide-sense stationary multivariate system observed through partial keyholes (sensor networks, distributed systems), and the core move — aggregate alignment-free marginal statistics, evaluate the global causal functional once — is the kind of framing that tends to spawn follow-ups (better completion priors, calibrated edge-level inference, GPU-scaled versions are all flagged as future work).

Where to spend your time: Section 2.3 and Fig. 1 for the mechanism (ten minutes and you could reimplement it), Section 3.6 for the three-way validation logic that makes the IBL claim credible, and Section 3.7 with Fig. 6 for the theta hierarchy — the part that turns a methods paper into a neuroscience result.