Non-stationary dynamics of interspike intervals in neuronal populations

TL;DR

The paper addresses how interspike interval statistics evolve in neuronal populations under nonstationary input by formulating a joint diffusion process in the membrane potential and time since last spike, leading to a two-dimensional Fokker-Planck equation for the population density. It develops a one-dimensional hierarchy of ISI moments via a τ-transform, and derives explicit stationary solutions and relaxation dynamics, as well as a first-order linear response to weak input modulation. The approach is validated against large-scale spiking-network simulations, showing accurate capture of ISI dynamics even far from stationarity and during limit-cycle oscillations, and provides a practical, computationally efficient tool through the moment hierarchy. This framework extends renewal-based theories, connects to age-structured dynamics, and offers a general method for analyzing threshold-crossing statistics in driven excitable systems with self-consistent inputs.

Abstract

We study the joint dynamics of membrane potential and time since the last spike in a population of integrate-and-fire neurons using a population density framework. This leads to a two-dimensional Fokker-Planck equation that captures the evolution of the full neuronal state, along with a one-dimensional hierarchy of equations for the moments of the inter-spike interval (ISI). The formalism allows us to characterize the time-dependent ISI distribution, even when the population is far from stationarity, such as under time-varying external input or during network oscillations. By performing a perturbative expansion around the stationary state, we also derive an analytic expression for the linear response of the ISI distribution to weak input modulations.

Figures (4)

• Figure 1: Relaxation to the stationary state in drift-dominated regime for a network of uncoupled leaky integrate-and-fire neurons (LIF, f(v) = -v): $\mu_0=21 \, \text{mV}$, $D_0=14.15 \, \text{mV}^2$, $H=0 \, \text{mV}$, $\theta=20 \, \text{mV}$. (a) Top: time-dependent ISI distribution, comparing Eq. \ref{['eq:relax-H']} with spiking neural network simulation (SNN $N=10^4$ neurons). Bottom: population density $q(v,\tau,t)$. (b) Temporal dynamics of mean ISI, comparison between theory and SNN simulation. (c) Temporal trajectory of mean ISI and coefficient of variation. Comparison between theory and SNN simulation.
• Figure 2: Linear response and transfer functions of first ISI moments. (a) Transfer function for the first two moments. Top: drift-dominated regime (same parameters as Fig. \ref{['fig:relax']}) and noise-dominated regime ($\mu = 16.314 \, \text{mV}/\tau_m$, $D_0=14.25 \, \text{mV}^2/\tau_m$, $H=10 \, \text{mV}$, $\theta=20 \, \text{mV}$). comparison between theory (solid lines) and moment hierarchy numerical integration Eq. \ref{['eq:moments']}. Solid grey line represents the firing rate transfer function $\mathrm{H}^\nu$. (b) Linear response theory for the ISI distribution in drift-dominated regime, $\omega/2\pi = 0.4\tau_m$. Left: drift modulation, right: noise modulation. Comparison between theory and numerical integration of Eq. \ref{['eq:isi-fp']}.
• Figure 3: Synchronous (limit cycle) irregular state in an excitatory population of LIF neurons. (a) ISI distribution $\rho_0(\tau) = \langle \nu_q(\tau,t)\rangle_t/\langle \nu(t)\rangle_t$. Comparison between population dynamics, numerical integration of Eq. \ref{['eq:isi-fp']} and SNN simulation ($N=5\cdot10^4$). (b) Phase-dependent ISI distribution. (c) Joint trajectory of firing rate, ISI mean (left) and coefficient of variation (right). Comparison between population dynamics \ref{['eq:isi-fp']}, moment hierarchy from Eq. \ref{['eq:moments']}) and SNN simulation ($N=5\cdot10^4$). Parameters: $J=0.1 \, \text{mV}$, $K=1000$, $\mu_\text{ext} = 21 \, \text{mV}/\tau_m$, $D_\text{ext}=1.33 \, \text{mV}/\tau_m$, $H=0 \, \text{mV}$, $\theta=20\, \text{mV}$, $\delta = 0.05 \, \tau_m$, $\tau_\delta = 0.1 \tau_m$.
• Figure 4: Synchronous (limit cycle) irregular in an inhibitory population of LIF neurons. (a) ISI distribution $\rho_0(\tau) = \langle \nu_q(\tau,t)\rangle_t/\langle \nu(t)\rangle_t$. Comparison between population dynamics, numerical integration of Eq. \ref{['eq:isi-fp']}) and SNN simulation ($N=5\cdot10^4$). (b) Time trajectory of mean ISI. Comparison between population dynamics \ref{['eq:isi-fp']}, moment hierarchy \ref{['eq:moments']} and SNN simulation ($N=5\cdot10^4$). (c) Same as (b) but tracking squared coefficient of variation. Parameters: $J=-0.1 \, \text{mV}$, $K=1000$, $D_\text{ext} = 8 \, \text{mV}$, $\delta = 0.4 \, \tau_m$, $\tau_\delta = 0\, \tau_m$. Other parameters as in Fig. \ref{['fig:3']}.