Travelling waves modulated by subthreshold oscillations in networks of integrate-and-fire neurons

Basic Information

• Paper ID: 2511.05232
• Title: Travelling waves modulated by subthreshold oscillations in networks of integrate-and-fire neurons
• Authors: Henry D. J. Kerr, Peter Ashwin, Kyle C. A. Wedgwood (University of Exeter)
• Classification: q-bio.NC (Quantitative Biology - Neurons and Cognition)
• Publication Date: November 7, 2025 (arXiv preprint)
• Paper Link: https://arxiv.org/abs/2511.05232

Abstract

This study investigates travelling wave phenomena in neural networks, specifically examining the effects of incorporating slow voltage-gated ion channels as linear adaptation variables into leaky integrate-and-fire (LIF) neural networks with lateral inhibition. The research reveals: (1) increased ion channel strength enhances wave speed through interaction with lateral inhibition; (2) ion channel-induced subthreshold oscillations structure the interspike intervals of weakly-coupled travelling wave solutions around the neuron's intrinsic oscillatory period; (3) subthreshold oscillations generate codimension-2 grazing bifurcations, producing fold bifurcations on wave solution branches and introducing multiple intermediate structures in solution space. The study also develops novel root-finding algorithms to handle uncertainties in the existence of spiking solutions.

Research Background and Motivation

Research Questions

This paper investigates the propagation mechanisms of travelling waves in neural networks, particularly how subthreshold oscillations at the single-neuron level influence collective behavior at the network level.

Problem Significance

1. Physiological Relevance: Travelling waves are ubiquitous in sensory, motor, and cognitive processes in the brain and are considered fundamental to brain computation.
2. Cross-Scale Dynamics: Subthreshold oscillations and resonant responses in local dynamics significantly affect travelling wave behavior. For example, in cochlear auditory processing, resonant frequencies directly relate to wave propagation speed.
3. Theoretical Value: Understanding how single-neuron dynamics shape network-level collective behavior is a core question in neuroscience.

Limitations of Existing Approaches

1. Simplified Models: Amari-type rate models oversimplify local dynamics, neglecting single-neuron and small-circuit processing (e.g., resonant responses).
2. Complex Models: Hodgkin-Huxley-type models, while detailed, are analytically intractable.
3. Lack of Subthreshold Dynamics: Previous IF network studies (e.g., reference 12) primarily focused on one-dimensional local dynamics, overlooking subthreshold oscillations.

Research Motivation

The authors seek a balance between the simplicity of IF models and biological realism by introducing linear adaptation variables (representing HCN or Kv1 ion channels) to capture subthreshold oscillations while maintaining analytical tractability.

Core Contributions

1. Model Extension: Extends previous LIF network models to include two-dimensional local dynamics with linear ion channel variables capable of generating subthreshold oscillations.
2. Analytical Construction: Establishes semi-explicit construction methods for travelling wave solutions and linear stability analysis frameworks (equations 34 and 46).
3. Bifurcation Structure Discovery:
   • Identifies codimension-2 double grazing bifurcation points (Type III grazing bifurcation)
   • Discovers complex bifurcation structures driven by subthreshold oscillations
   • Reveals distinctions between "atomic waves" and "composite waves"
4. Wave Speed Modulation Mechanism: Clarifies how ion channel parameters R (response rate) and D (decay rate) influence wave speed through interaction with lateral inhibition.
5. Locking Phenomena: Discovers that interspike intervals of weakly-coupled biphasic waves lock to integer multiples of the neuron's intrinsic period.
6. Efficient Numerical Algorithm: Develops event-driven GPU-accelerated simulation algorithms using improved Newton-Raphson methods to precisely capture spike times.

Detailed Methods

Task Definition

Input: N uniformly distributed LIF neurons on a ring domain with Mexican hat connectivity kernel

Output: Existence, speed, stability of travelling wave solutions and their bifurcation structure as parameters vary

Constraints: Neurons coupled through distance-dependent all-to-all connections, following leaky integrate-and-fire dynamics

Model Architecture

1. Single Neuron Dynamics

For the n-th neuron (n=1,...,N), the model is described by three ordinary differential equations:

$\frac{dv_n}{dt} = I - v_n - u_n + s_n - (v_{th} - v_r)\sum_{k\in\mathbb{Z}}\delta(t-t_{n,k})$ (1)

$\frac{du_n}{dt} = Rv_n - Du_n$ (2)

$\frac{ds_n}{dt} = -\beta s_n + \beta f^{in}_n(t)$ (3)

Where:

• $v_n$: membrane potential (spikes and resets to $v_r=0$ when reaching threshold $v_{th}=1$)
• $u_n$: ion channel current (linearized adaptation variable)
• $s_n$: synaptic buffer variable
• $R\geq 0$: response rate of ion channel to voltage changes
• $D>0$: decay rate of ion channel
• $\beta>0$: synaptic timescale parameter
• $I$: external input current (controls excitability)

Key Innovation: The variable $u_n$ models voltage-gated ion channels through parameters R and D, capable of generating subthreshold oscillations. The system matrix eigenvalues are $\lambda_{1,2} = -p\pm q$, where: $p = \frac{1}{2}(D+1), \quad q = \frac{1}{2}\sqrt{(D-1)^2 - 4R}$

When $4R > (D-1)^2$, q is imaginary, and the system exhibits damped oscillations with intrinsic frequency $|q|/2\pi$.

2. Connectivity Structure

Mexican hat connectivity kernel: $w(d) = \frac{A}{a\sqrt{2\pi}}e^{-\frac{d^2}{2a^2}} - \frac{B}{b\sqrt{2\pi}}e^{-\frac{d^2}{2b^2}}$ (7)

Where $A=B$ ensures excitatory-inhibitory balance, and $b>a$ produces short-range excitation, intermediate-range inhibition, and long-range no interaction.

Synaptic input: $f^{in}_n(t) = \Delta x\sum_{k\in\mathbb{Z}}\sum_{m\neq n}w(|x_n-x_m|)\delta(t-t_{m,k})$ (6)

3. Continuum Limit and Travelling Wave Solution Construction

Taking the $N\to\infty$ continuum limit, introduce spike time function $t_k(x)$ (k-th spike) and spike position function $X_k(t)$.

Travelling wave solution form: $t_j(x) = \tau_j + \frac{x}{c}, \quad j=1,2,...,m$ (26)

Where c is wave speed and $\tau_j$ are time offsets (setting $\tau_1=0$).

In the co-moving coordinate system $\xi = t - x/c$, the system becomes time-invariant, and travelling wave solutions satisfy: $\mathbf{v}(\xi) = Ie^{\xi M}\int_{-\infty}^{\xi}e^{-\zeta M}\begin{pmatrix}1\\0\end{pmatrix}d\zeta + \sum_{j=1}^m\beta e^{\xi M}\int_{-\infty}^{\xi}e^{-\zeta M}\begin{pmatrix}1\\0\end{pmatrix}e^{-\beta\zeta}\int_{-\infty}^{\zeta}e^{\beta r}w(c(r-\tau_j))c\,dr\,d\zeta$$- (v_{th}-v_r)\sum_{j=1}^m e^{(\xi-\tau_j)M}\begin{pmatrix}1\\0\end{pmatrix}\Theta(\xi-\tau_j)$ (34)

The m spike events provide m conditions $v(\tau_j^-) = v_{th}$, solved using Newton-Raphson methods for $(c, \tau_2,...,\tau_m)$.

4. Stability Analysis

Introduce perturbations $\tilde{X}_j(t) = c(t-\tau_j) + \epsilon\phi_j(t)$, where $\phi_j(t) = \text{Re}(\Phi_j e^{\lambda t})$.

After linearization, the characteristic equation becomes: $\det(F(\lambda) - G) = 0$ (46)

Where $F(\lambda)$ is an $m\times m$ matrix and $G$ is a diagonal matrix with elements defined by equation (43). The wave is stable if all $\lambda$ (except $\lambda=0$ from translational invariance) have negative real parts.

Technical Innovations

1. Parameter Choice: Uses $(R,D)$ parameterization rather than traditional conductance parameters, facilitating visualization and comparison while clearly connecting parameter space to intrinsic oscillations and resonant responses.
2. Semi-Explicit Solutions: Exploits system linearity to construct explicit solutions between spike events (equations 12-13), avoiding cumulative numerical integration errors.
3. Improved Newton-Raphson Algorithm:
   • Constructs upper bounds $m_n \geq \sup\{v'(t): t\in[t_n, t_{n+1}]\}$ to prevent overshoot
   • Combines two bounds achieving fast convergence away from roots and quadratic convergence near roots
   • Detects solution non-existence (when $t_n > T$ or $M_n \leq 0$)
4. Event-Driven Simulation: Directly jumps from one spike event to the next, utilizing GPU parallelization to compute next spike times for each neuron.

Experimental Setup

Model Parameters

Baseline parameter set (Table 1):

• Ion channel decay rate: $D = 1$
• Synaptic response rate: $\beta = 6$
• Resting voltage: $v_{rest} = 0.9$
• Connection strength: $A = B = 2$
• Connection range: $a = 1, b = 2$
• Number of neurons: $N = 2000$
• Domain length: $2L = 20$

Varied Parameters

Two main parameter spaces explored:

1. (R, D) space: Ion channel parameters, focusing on oscillatory region $4R > (D-1)^2$
2. ($\beta$, R) space: Synaptic timescale versus ion channel response rate

Numerical Methods

1. Pseudo-arclength continuation: Tracks solution branches as parameters vary
2. Bifurcation detection: Identifies grazing, fold, and Hopf bifurcations
3. Slow parameter variation simulations: $R(t) = R_0 + \delta_R \min\{t, t_{fin}\}$ to study dynamics crossing bifurcations

Evaluation Metrics

• Wave speed c
• Interspike times $\tau_2, \tau_3,...$
• Stability spectrum (real parts of eigenvalues $\lambda$)
• Bifurcation point locations

Experimental Results

Main Results

1. Wave Speed under R Variation (Figure 5, D=1, β=6)

Uniphasic Waves:

• Wave speed increases monotonically with R (from c≈1.5 at R=0 to c≈3.5 at R=4)
• Grazing bifurcation occurs at R≈1.9, rendering solutions inadmissible

Biphasic Atomic Waves (Slow Branch):

• Wave speed approximately 0.5 units lower than uniphasic waves
• Single s peak between v peaks with monotonic v rise
• Grazing bifurcation terminates at R≈1.87
• Speed difference from uniphasic waves essentially constant (validating transient response intuition from Section 2.5)

Biphasic Composite Waves (Fast Branch):

• Wave speed close to uniphasic waves
• Longer interspike time $\tau_2$ with two s peaks
• Fold bifurcation at R≈2.8 producing stable-unstable pairs
• Unstable waves slowly converge to stable solutions in simulations (delayed bifurcation in Example 2)

2. Locking Phenomena in Weakly-Coupled Waves (Figure 7)

Multiple biphasic wave branches discovered with $\tau_2$ values near integer multiples of intrinsic period $2\pi/|q|$:

• D=0.85: Odd multiples stable for R<5; even multiples for R>5
• D=0.88: Transition state with branch reconfiguration at R≈5, fold connections observed
• D=1: Even multiples stable for R<5; odd multiples for R>5

Stability exhibits alternating odd-even patterns, indicating weakly-coupled biphasic waves are essentially "locked" combinations of two uniphasic waves.

3. Two-Parameter Bifurcation Diagram (Figure 8, R-D plane)

Key structures identified:

• Double Grazing Points: Grazing curves of uniphasic and biphasic slow branches converge at specific (R,D) values, corresponding to two local maxima simultaneously touching threshold (codimension-2 bifurcation)
• Fold Bifurcation Curves: Fast biphasic branch fold curves partition parameter space into solution existence/non-existence regions
• Oscillation Boundary: $4R = (D-1)^2$ curve marks the boundary between triangular/hyperbolic solutions

General pattern: Strong adaptation variables (high R, low D) eliminate these wave solutions, while timescale changes have smaller effects.

4. β Parameter Variation (Figure 10, R=2.5, 2.6, 2.7)

• Admissible solutions exist only at higher β values, with wave speed increasing with β
• Increasing R from 2.6 to 2.7 produces two fold bifurcations in fast biphasic branch, creating intermediate β gaps (solution non-existence)
• Compared to R=0 previous results: R>0 introduces fold bifurcations, while R=0 primarily exhibits Hopf bifurcations

5. Bifurcation Crossing Simulations (Figures 9 and 10)

Grazing Bifurcation (Figures 9.1 and 9.3):

• Slow biphasic waves: entire wave immediately transforms to single bump
• Fast biphasic waves: second wave component forms bump, first component continues propagating until suppressed, forming second bump

Fold Bifurcation (Figures 9.2 and 10):

• Exhibits delayed bifurcation characteristics, waves continue propagating well beyond bifurcation point
• Eventually transform to bumps, which may be unstable and split into multiple bumps
• Wave reflection phenomena and transient bistability with tonic firing observed

Ablation Studies

While not explicitly labeled as "ablation studies," the paper conducts substantive component analysis through systematic single-parameter variation:

1. R=0 Baseline: Comparison with reference 12 results validates recovery of original behavior when R=0
2. Fixed D, Varying R: Isolates ion channel response rate effects
3. Fixed R, Varying D: Isolates timescale effects
4. Fixed R, Varying β: Isolates synaptic timescale effects

Case Studies

Example 1 (Figure 5, R=3.5): Stable fast biphasic wave

• Contours show two clearly separated s peaks
• Stable wave propagation in simulations with regular spike time diagonal lines

Examples 4-6 (Figure 5, Unstable Waves):

• Examples 4 and 5: Immediately lose one peak, oscillate then converge to uniphasic wave
• Example 6: Near grazing point, instability produces additional spike events propagating in both directions, eventually network quiescence

Slow Parameter Variation (Figure 9.2):

• Wave persists approximately 60 time units before transforming across fold bifurcation
• Bump formation exhibits chevron pattern

Experimental Findings

1. Wave Speed Modulation Mechanism: R increase → faster neuron response → increased wave speed (consistent with transient response analysis)
2. Peak Number and Speed: Atomic m-peak waves slower than uniphasic waves, with speed decreasing as m increases (Figure 6)
3. Stability Transitions: Higher m-value waves lose stability at lower R values (m=7 stable region disappears)
4. Locking-Unlocking Transitions: Locking patterns of weakly-coupled waves reorganize at specific R values
5. Bifurcation Cascades: Atomic waves undergo Hopf→grazing bifurcation sequences as R increases
6. Double Grazing Organization: Acts as parameter space organizational center, separating regions of different grazing mechanisms

Related Work

Neural Field Models

1. Amari (1977): Pioneering rate model establishing theoretical foundations for travelling waves and bumps in lateral inhibition networks
2. Ermentrout, Bressloff, Coombes, et al.: Developed neural field theory, analyzing effects of synaptic delays, recurrent inhibition, anisotropy on wave propagation

Travelling Waves in Spiking Networks

1. Laing & Chow (2001): First construction of bump attractors in LIF networks
2. Avitabile, Davis & Wedgwood (2023): Systematic study of bumps and travelling waves in R=0 LIF networks, discovering saddle waves and spatiotemporal chaos
3. Bressloff (2000): Established stability analysis framework for travelling waves in excitatory IF networks

Subthreshold Oscillations and Resonance

1. Richardson, Brunel & Hakim (2003); Rotstein & Nadim (2014): Quantified resonance mechanisms in IF neurons
2. Stark et al. (2022): Demonstrated network resonance can independently arise at multiple neural organizational levels
3. Bell (2004, 2012); Nankali et al. (2022): Relationships between travelling waves and resonance in cochlea

Advantages of This Work

• Compared to Amari models: Retains single-neuron dynamics details, revealing subthreshold oscillation effects on network behavior
• Compared to Hodgkin-Huxley models: Maintains analytical tractability, establishes semi-explicit solutions
• Compared to Avitabile et al. (2023): Extends to two-dimensional local dynamics, discovering new bifurcation structures (double grazing, locking phenomena)
• Compared to resonance studies: First systematic analysis of subthreshold oscillations in spiking network travelling waves

Conclusions and Discussion

Main Conclusions

1. Slow ion channels increase wave speed through interaction with lateral inhibition, with effects monotonically increasing with R
2. Subthreshold oscillations lock interspike intervals of weakly-coupled waves to integer multiples of intrinsic period, forming phenomena analogous to convection cell tail binding
3. Double grazing bifurcation serves as parameter space organizational center, producing codimension-2 Type III grazing (Class A discontinuous systems)
4. Atomic and composite waves exhibit different bifurcation behavior: former primarily undergoes grazing, latter involves fold bifurcations
5. Bifurcation crossing dynamics depend on bifurcation type: grazing causes rapid transitions, folds produce delayed bifurcations
6. Wave-to-bump transition mechanisms are diverse: can form single or multiple bumps, potentially accompanied by wave reflection

Limitations

1. Simplified Connectivity: Mexican hat kernel merges excitation and inhibition, potentially missing important dynamics from separated E-I populations (e.g., slow waves with inhibitory cells firing before excitatory cells)
2. Dale's Principle: Single neurons producing both excitation and inhibition violates biological principles
3. Linear Ion Channels: Actual ion channels exhibit nonlinear conductances, limiting precise modeling of specific channel types
4. One-Dimensional Domain: Unexplored are planar wave instabilities, spiral waves, gliders in higher dimensions
5. Deterministic Model: Ignores stochasticity, cannot explain waves with partial neuron participation
6. Homogeneous Networks: Neglects heterogeneity and anisotropy

Future Directions

1. Extension to Three-Dimensional Local Dynamics: Introduce ultra-slow timescales producing bursting, potentially using piecewise-linear models maintaining tractability
2. Separated E-I Populations: Respect Dale's principle, study bump interactions between different populations
3. Dynamic Input Conductance Models: Integrate fast, slow, ultra-slow components, connecting specific ion channel expression to network behavior
4. High-Dimensional Extensions: Analyze wave front instabilities in two and three dimensions, study gliders and localized travelling wave structures
5. Stochasticity and Heterogeneity: Use Hawkes processes describing spike propensities, investigate lurching waves and complex spatiotemporal patterns
6. Cross-Scale Oscillations: Study mechanisms of up-down state transitions in slow cortical waves (e.g., sleep slow oscillations)
7. Neural Fields on Surfaces: Extend to arbitrary surfaces (e.g., cortical surfaces, organoids) with geometric structure

In-Depth Evaluation

Strengths

1. Theoretical Depth:
   • Establishes complete analytical framework (construction + stability), extending Bressloff (2000) methods
   • Discovers codimension-2 double grazing bifurcation as novel organizational structure
   • Reveals physical mechanisms of locking phenomena (integer multiples of intrinsic period)
2. Methodological Innovation:
   • Improved Newton-Raphson algorithm elegantly handles spike time uncertainties
   • Event-driven simulation avoids inefficiency and errors of traditional time-stepping
   • GPU parallelization fully exploits problem's "embarrassingly parallel" nature
3. Systematicity:
   • Comprehensive exploration of (R,D) and (β,R) parameter spaces
   • Tracks atomic wave branches up to 10 peaks
   • Combines analytical, numerical continuation, and direct simulation methods for mutual validation
4. Physical Intuition:
   • Section 2.5 transient response analysis provides clear intuition for main results
   • Atomic vs. composite wave conceptualization aids understanding of branch essences
5. Reproducibility:
   • Open-source code (GitHub)
   • Complete parameter tables (Table 1)
   • Detailed algorithm pseudocode (Appendix B)

Weaknesses

1. Biological Realism:
   • Linear ion channels are crude approximations
   • Mexican hat connectivity violates Dale's principle
   • Lacks noise and heterogeneity
2. Analytical Limitations:
   • High-peak waves (m>10) insufficiently explored
   • Three-dimensional local dynamics (e.g., bursting) exceed current framework
   • Bump stability not systematically analyzed (only observed through simulation)
3. Missing Experimental Validation:
   • No comparison with experimental data
   • Biological justification of parameter choices insufficiently argued
   • Predictive testability unclear
4. Incomplete Phenomenon Explanations:
   • Criteria for wave reflection vs. bump formation not explicit
   • Mechanisms of transient bistability with tonic firing not deeply analyzed
   • Complete theory of locking pattern reorganization (Figure 7) lacking
5. Writing Details:
   • Some symbols undefined at first appearance (e.g., $\Theta$ in equation 34)
   • Figure 5's six examples have extremely high information density, requiring repeated reference to text

Impact

Contributions to the Field:

• Fills theoretical gap between rate models and detailed spiking models
• Provides quantitative framework for understanding subthreshold dynamics effects at network level
• Double grazing bifurcation discovery may inspire research in other excitable systems

Practical Value:

• Numerical algorithms have practical application value for large-scale IF network simulations
• Provides guidance for designing neuromorphic hardware with specific wave speeds
• May help understand pathological/physiological phenomena (epileptic waves, sleep slow waves)

Theoretical Impact:

• Demonstrates IF networks can exhibit behavior qualitatively different from continuous models (e.g., spatiotemporal chaos of bumps)
• Locking phenomena and fluid mechanics convection cell analogy may promote interdisciplinary exchange
• Detailed characterization of bifurcation cascades and delayed bifurcations enriches non-smooth dynamical systems theory

Reproducibility: Excellent

• Open-source code
• Detailed algorithms
• Complete parameters
• Likely to become benchmark implementation for this direction

Applicable Scenarios

1. Theoretical Neuroscience: Studying contributions of specific ion channels (HCN, Kv1) to network oscillations
2. Computational Neuroscience: Large-scale network simulations balancing efficiency and biological realism
3. Neuromorphic Engineering: Designing spiking computational architectures utilizing wave propagation for information processing
4. Dynamical Systems Theory: Case studies for excitable media, non-smooth systems, delayed bifurcations
5. Medical Applications: Understanding potential mechanisms of pathological neural waves (epileptic seizure propagation)

Inapplicable Scenarios:

• Research requiring precise matching of specific experimental data (model too simplified)
• Studying synaptic plasticity, learning, and other long-timescale phenomena (static connectivity model)
• Pharmacological research requiring detailed ion channel dynamics

Key References

1. Amari (1977): Foundational work in neural field theory
2. Bressloff (2000): Pioneering IF network travelling wave stability analysis methods
3. Avitabile, Davis & Wedgwood (2023): Direct baseline study (R=0) extended by this work
4. Richardson, Brunel & Hakim (2003): Quantitative theory of subthreshold resonance
5. Kowalczyk et al. (2006): Classification framework for double grazing bifurcations
6. Laing & Chow (2001): First construction of bump attractors in spiking networks
7. Ermentrout (1998): Classical analysis of synaptically-generated travelling waves

Overall Assessment: This is a high-quality theoretical neuroscience paper, rigorous and innovative in methodology, insightful in results. By introducing subthreshold oscillations as a biologically relevant feature, the authors discover rich bifurcation structures and locking phenomena, significantly advancing our understanding of collective behavior in spiking neural networks. Open-source efficient numerical tools further enhance practical value. Main limitations lie in model simplifications and lack of experimental validation, but as theoretical exploration, this work achieves high standards and is expected to exert sustained impact on neural dynamics and excitable media research.