Statistical Inference and Stability Boundaries of Multi-cellular Interaction Hypergraphs from Asynchronous Event Streams

Abstract: We introduce the Hyperedge-triggered Hawkes (HTH) process for inferring higher-order interaction structure in multi-cellular systems from asynchronous event-time data. Beyond standard pairwise excitation, the HTH intensity includes a term activated by the simultaneous co-firing of a cell group within a temporal window. We derive a closed-form Expectation-Maximisation algorithm whose key ingredient is a piecewise compensator that eliminates the systematic bias present in the naive integral formulation. A CP tensor decomposition reduces the hyperedge parameter count from O(N^K) to O(NR). Across eleven synthetic experiments the framework achieves pairwise recovery error below 5%, while revealing a systematic -22% bias on hyperedge weights that is non-monotonic in the kernel decay rate, ruling out a simple temporal-overlap explanation and motivating adaptive kernel methods. On multi-electrode recordings of mouse retinal ganglion cells, the model yields a +20.6 nat likelihood gain over the pairwise baseline, providing suggestive but not decisive evidence for higher-order interactions. Code and all experiments are publicly available at https://github.com/Hanii0210/hypergraph-hawkes.

• The paper introduces a hyperedge-triggered Hawkes process with pattern-completion anchors, a corrected piecewise compensator, CP tensor parameterization, and closed-form EM updates for scalable higher-order interaction inference.
• Synthetic experiments recover baseline and pairwise parameters within 5%, identify stability near a spectral radius of 1, and show reliable detection of true hyperedges, but hyperedge weights remain systematically underestimated by 22%.
• The paper finds suggestive higher-order structure in seven-neuron retinal recordings, with a likelihood gain of 20.6 nats, although BIC favors the simpler pairwise model and highlights risks from short recordings, heterogeneous firing rates, and stimulus-driven correlations.

The Hyperedge-triggered Hawkes (HTH) process introduced by Xu extends the classical mutually-exciting point process framework to represent genuinely higher-order interactions among cells, where the near-simultaneous co-firing of a group triggers downstream activity beyond the sum of pairwise contributions. The paper contributes a mathematically specified intensity model, a closed-form Expectation–Maximisation (EM) algorithm with a corrected piecewise compensator, a CP tensor parameterisation for tractability, and an extensive empirical evaluation spanning eleven synthetic benchmarks and one real neural recording dataset. Notably, the work is candid about its own failure modes: it documents a systematic $-22\%$ bias on inferred hyperedge weights and reports that on real data the likelihood favours higher-order structure while BIC does not.

Model formulation

The conditional intensity of node $n$ comprises three additive terms: a baseline rate $\mu_n$, standard pairwise Hawkes excitation $\alpha_{n_j \to n}\,\phi(t - t_j)$ with exponential kernel $\phi(\tau) = e^{-\beta \tau}$, and a hyperedge term activated by pattern completion. A pattern completion for hyperedge $e$ occurs when every member has fired within a temporal window $\Delta$; the anchor time is the most recent such completion before $t$, selected via a $\max$ operator. This "pattern-completion anchor" semantics means only the latest completion contributes to the current intensity, which prevents double-counting across overlapping completions and keeps the likelihood tractable — at the cost of discarding information from earlier completions within the same window.

Hyperedge weights are parameterised through a rank-$R$ CP decomposition over a non-negative factor matrix $n$0, reducing the free parameter count from $n$1 to $n$2 without restricting which interaction patterns are representable.

A key methodological contribution is the piecewise compensator. Because each anchor remains active only until the next pattern completion, integrating each anchor's kernel contribution from activation time all the way to $n$3 systematically overcounts and biases $n$4 downward. The correct compensator sums truncated kernel integrals between consecutive completion times:

$n$5

The author reports catching this bug during development and validating the fix experimentally; the residual bias discussed below persists even after the correction.

Inference

The E-step decomposes each event's responsibility across background firing, pairwise excitation from a specific prior event, and triggering by a specific hyperedge, with responsibilities normalised by the total intensity. All three parameter types admit closed-form M-step updates: the baseline is a responsibility-weighted event count divided by $n$6; pairwise weights use a standard truncated-kernel denominator; hyperedge weights divide attributed mass by the piecewise compensator plus an $n$7 penalty enforcing sparsity.

Because enumerating all $n$8 candidate hyperedges is combinatorially prohibitive, inference proceeds via a two-stage pipeline: a pairwise-only Hawkes fit flags significantly mutually exciting node pairs, candidate hyperedges are formed as cliques in this significance graph, and full HTH inference with $n$9 regularisation retains only surviving edges. This design implicitly assumes that genuine higher-order interactions leave a detectable pairwise footprint — a plausible but untested assumption that could cause the filter to miss hyperedges whose members are only weakly pairwise coupled.

Synthetic validation

All synthetic experiments use Ogata thinning simulators with known ground truth. Across 25 independent datasets on a canonical 3-node system, baseline rates and pairwise weights recover with relative error below 5%, while the hyperedge weight exhibits a systematic downward bias of $\mu_n$0 (mean inferred 0.314 vs. true 0.400) with coefficient of variation 0.36. The bias is consistent across seeds, indicating a structural rather than stochastic origin. This is the paper's most consequential negative result: the estimator is reliable for pairwise structure but materially underestimates group-level excitation strength.

Supporting experiments establish the framework's other properties. An $\mu_n$1 regularisation path retains the true hyperedge across two orders of magnitude of penalty while eliminating five decoys, with AIC and BIC agreeing on $\mu_n$2. Twenty EM runs from random initialisations converge to the same log-likelihood (std = 0.004 nats), consistent with a unique optimum. As hyperedge strength grows, the spectral radius $\mu_n$3 of the combined interaction matrix crosses unity at approximately 1.24 times the theoretically predicted critical strength, marking the transition from stable dynamics to unbounded cascades; the 24% offset aligns closely with the $\mu_n$4 recovery bias, suggesting a shared mechanistic cause. An independent copula-based analysis confirms that hyperedge interactions leave a detectable fingerprint in upper-tail dependence ($\mu_n$5 against the pairwise-only null). On a 4-node system with a true 3-node hyperedge, recovery error is 9% with all decoys suppressed below 0.01, confirming generalisation beyond order two.

Two falsifiability-oriented results deserve emphasis. In a symmetric test, data generated with a true hyperedge yield a likelihood gain of $\mu_n$6 (BIC difference +10.97), decisively favouring HTH, whereas data generated without one yield $\mu_n$7 (BIC difference −6.27): the method finds hyperedges when they exist and does not invent them when they do not. Sweeping the window parameter $\mu_n$8 shows a sharply peaked log-likelihood at the true value, where the hyperedge weight error drops to 0.7%, establishing identifiability of $\mu_n$9 by simple grid search. Wall-clock time per EM iteration scales as $\alpha_{n_j \to n}\,\phi(t - t_j)$0, matching the theoretical $\alpha_{n_j \to n}\,\phi(t - t_j)$1 bound, with per-event-pair cost around 2.2 μs.

Bias ablation

To probe the bias mechanism, the author sweeps the kernel decay rate $\alpha_{n_j \to n}\,\phi(t - t_j)$2 over five values with all else fixed. If the bias arose purely from temporal overlap between pairwise and hyperedge contributions, it should grow monotonically as decay slows. It does not: the bias is −18.9% at $\alpha_{n_j \to n}\,\phi(t - t_j)$3, nearly zero at intermediate rates (−2.9% at $\alpha_{n_j \to n}\,\phi(t - t_j)$4, −0.5% at $\alpha_{n_j \to n}\,\phi(t - t_j)$5), and flips positive (+12.3%) at fast decay, where variance also explodes (std/mean > 1), consistent with overfitting under a too-narrow kernel window. This non-monotonicity rules out the simple temporal-overlap explanation and implicates an interaction between component coupling and kernel misspecification. The practical implication is that any deployed version of HTH requires adaptive or per-node kernel bandwidth selection; a single global $\alpha_{n_j \to n}\,\phi(t - t_j)$6 yields biased hyperedge estimates whose sign depends on the timescale regime.

Real-data application

On multi-electrode recordings of 7 mouse retinal ganglion cells (CRCNS ret-1), using a 40-second window with 3,759 spikes and firing rates from 3 Hz to 29 Hz, the full HTH model (47 parameters, 6 candidate hyperedges) reaches a log-likelihood of 6780.9 versus 6760.2 for the pairwise-only model: a gain of $\alpha_{n_j \to n}\,\phi(t - t_j)$7 nats, giving a likelihood-ratio statistic of 41.2 against a $\alpha_{n_j \to n}\,\phi(t - t_j)$8 critical value of 12.59 for 6 degrees of freedom. Five of six candidates survive with weights between 0.06 and 0.35, and neuron 0 emerges as a hub receiving strong pairwise excitation ($\alpha_{n_j \to n}\,\phi(t - t_j)$9) from neurons 3, 4, and 6, consistent with known retinal architecture.

The BIC verdict, however, reverses: BIC difference of −8.1 favours the simpler pairwise model. The paper states plainly that "the likelihood says yes; BIC says not yet," characterising the evidence as suggestive rather than decisive. Three factors widen the gap: ten-fold variation in firing rates violates the assumption of comparable baselines; the binary white-noise stimulus can drive correlated responses that mimic hyperedge interactions without reflecting synaptic group coupling; and 40 seconds of data is short. One candidate, $\phi(\tau) = e^{-\beta \tau}$0, produces an anomalous weight of 3.73 despite neuron 5 having only 125 spikes — a sparse-data artifact in which low spike counts shrink the compensator integral and inflate the estimate. This artifact is a concrete caution for practitioners applying the method to heterogeneous-rate populations.

Limitations and open questions

The paper is explicit about its boundaries. The systematic $\phi(\tau) = e^{-\beta \tau}$1 hyperedge bias lacks a closed-form correction; the proposed remedy, a simulation-calibrated correction factor, would bring estimates within Monte Carlo error of truth but requires joint treatment of component separation and kernel selection given the non-monotonic $\phi(\tau) = e^{-\beta \tau}$2-dependence. The single global exponential kernel cannot accommodate diverse neuronal temporal profiles (transient versus sustained cells), motivating per-node decay rates $\phi(\tau) = e^{-\beta \tau}$3. The two-stage candidate generation assumes pairwise-detectable precursors of higher-order structure. On the real data, two interpretations remain unresolved: either genuine hyperedge interactions exist but the recording is too short for BIC to overcome the complexity penalty, or the observed co-activation is stimulus-driven rather than network-driven — distinguishable only through stimulus-conditional modelling. Additional open items include single-use anchor semantics to capture refractory dynamics, and computational scaling to $\phi(\tau) = e^{-\beta \tau}$4–$\phi(\tau) = e^{-\beta \tau}$5 events via vectorised or GPU implementations.

Conclusion

This paper delivers a complete, falsifiable pipeline for inferring higher-order interaction structure from asynchronous event-time data: a hyperedge-triggered intensity with pattern-completion anchors, a closed-form EM algorithm with a piecewise compensator correction, CP tensor parameterisation, and clique-based candidate generation. The synthetic evidence supports accurate pairwise recovery, stable convergence, sharp model selection, identifiable window parameters, and a detectable stability boundary at $\phi(\tau) = e^{-\beta \tau}$6. Equally valuable are the documented shortcomings — a structural $\phi(\tau) = e^{-\beta \tau}$7 hyperedge bias with non-monotonic kernel dependence, and a real-data result in which likelihood and BIC disagree. By reporting both strengths and failures with equal precision, the work provides a solid foundation for the adaptive-kernel and bias-correction methods its own results show to be necessary.