Exploring a regular-spiking limit of autoregressive point-process models with an absolute refractory period

These notes explore a deterministic limit of Autoregressive Point-Process Generalized Linear Model (AR-PPGLM) of a spiking neuron. We consider a Poisson point process that is either strongly driven or strongly suppressed, and take a limit where spiking output becomes deterministic (and instantaneous rates diverge). This deterministic limit is of interest in clarifying how AR-PPGLM models relate to deterministic models of spiking.

[get as PDF]

Introduction

Point-process models are used in neuroscience to understand what drives neurons to spike. Typically, these take the form of a Poisson process with a time-varying firing rate $\lambda (t)$ that is predicted from experimentally measured variables and the history of spiking activity.

The linear-nonlinear-Poisson formulation of point-process models doesn't capture the physiological conductances and currents that cause spiking. For this reason, these models typically only emulate neuronal activity seen in the training data . These models also treat spiking as a random process, but spiking can be quite regular at times. Models fit to regular spiking data learn to strongly increase or suppress their firing rate to capture this regularity, which this can lead to diverging parameters or rates.

In Rule and Sanguinetti (2008) , we explored these models mathematically, and showed one way to calculate how their mean rates and correlations evolve in time. These approximations assume that the underlying process is stochastic, with small Gaussian fluctuations. Considering what happens as we approach the regular spiking limit is therefore a stress-test of these models, to how (or whether) they start to break down.

An autoregressive point-process with an absolute refractory period

We model neuronal spike trains as a sequence of impulses (Dirac $\delta$ distributions) in time:

$$y(t)=\sum _i\delta (t-t_i),$$

where $t_i$ are the spike times. Autoregressive point-process GLMs predict the instantaneous firing rate $\lambda (t)$ of a neuron as a linear-nonlinear function of spiking history $h(t,\tau )$ and external input $u(t)$. One way to write a model for the output spike train $y(t)$ of a neuron is:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\lambda (t)& =\varphi (\text{H}h(t)+u(t))\\ y(t)& \sim \mathrm{Poisson}(\lambda (t)\cdot dt)\\ h(t,\tau )& =y(t-\tau ).\end{array}\end{array}$$

Here, $\varphi (\cdot )$ is a rectifying nonlinearity that reflects that the firing rate must be positive, and $\text{H}$ is a spiking history filter. The function $\varphi$ can be any rectifying nonlinearity, although it is usually exponential, logistic, or rectified polynomial. Above, we have denoted the history filter as a linear operator, defined for a single neuron as:

$$\text{H}h(y)=\underset{ϵ\downarrow 0}{lim}{\int }_{\mathrm{\infty }}^{ϵ}\eta (\tau )h(t,\tau )d\tau ,$$

where $\eta (\tau )$ are the history filter weights. The limit $ϵ\downarrow 0$ simply means that we exclude the present time-point ($\tau =0$) from the spiking history. In single neurons, the spike history filter captures phenomena like resonance, bursting, or a relative refractory period. In population models, it also captures how different neurons excite or inhibit one another.

After a neuron spikes, there is typically a brief 5-10 ms time window in which it cannot spike again, called the absolute refractory period. The firing rate is zero during this time. Depending on the choice of nonlinearity $\varphi$, it can be hard for the linear weights to drive the rate to zero. This can be fixed by adding a gating term that goes to zero if the cell has spiked recently:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\lambda (t)& =g[h(t)]\cdot \varphi (\text{H}h(t)+u(t))\end{array}\end{array}$$

We consider the case of a unit that is constrained to fire a most once every unit interval. That is, we express time in units such that the absolute refractory period is $\mathrm{\Delta }t=1$. The gating function $g[h(t)]$ is therefore defined as:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}g[h(t,\tau )]& =1-\underset{ϵ\downarrow 0}{lim}{\int }_{ϵ}^{1}h(t,\tau )d\tau .\\ & =1-\text{R}h(t,\tau )\end{array}\end{array}$$

This sets the rate to zero if a spike has occurred in the past $\tau \in (0,1]$ time-points.

Solution for steady-state firing rate in a simplified scenario

First we look at a mean-field model of the mapping between the external drive and output firing rate of a Poisson neuron with an absolute refractory period. Let ${\lambda }_0$ be the firing-rate response to a constant input for a neuron without the absolute refractory period. We neglect spiking entirely and assume that all dynamics summarized by the average rate $\lambda (t)$ that is constant in time. If a steady-state is stable, we can find for the steady-state rate with absolute refractoriness $\lambda$ by solving the relationship:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\lambda & =[1-{\int }_{\tau \in (0,1]}h(\tau )]{\lambda }_0{e}^{\lambda ⟨\eta ⟩}\\ & =[1-\lambda ]{\lambda }_0{e}^{\lambda ⟨\eta ⟩}.\end{array}\end{array}$$

Where $⟨\eta ⟩$ is the average value of this history filter. There is no closed-form solution to the above equation. If we ignore the history effects for now (or assume that they integrate to zero), we can solve the simpler relationship:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\lambda & =[1-\lambda ]{\lambda }_0,\end{array}\end{array}$$

which solves to $$\begin{array}{}\text{(6)}& \begin{array}{r}\lambda =\frac{{\lambda }_0}{1+{\lambda }_0}.\end{array}\end{array}$$

In the case of the Poisson GLM, our nonlinearity is exponential $\varphi =\exp$, and this reduces to a logistic function $\sigma (x)=[1+e^{-x}{]}^{-1}$:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\lambda & =\frac{{\lambda }_0}{1+{\lambda }_0}\\ & ={[1+{{\lambda }_0}^{-1}]}^{-1}\\ & ={[1+e^{-\log ({\lambda }_0)}]}^{-1}\\ & =\sigma [\log ({\lambda }_0)]\end{array}\end{array}$$

In this simplified scenario, then a Poisson GLM with an absolute refractory period behaves much like a discrete-time Bernoulli GLM, with a logistic sigmoid nonlinearity and binary spiking events $y_i\in \{0,1\}$ (Figure 1a). From a practical standpoint, this confirms the intuition that a good way to fit spike train data is to bin the spike train at a time resolution equal to the absolute refractory period, and fit a Bernoulli GLM. However, in the next section we will show that the observation model is not identical to the discrete-time Bernoulli GLM (Figure 1b).

In this scenario, the firing rate approaches $\lambda =1$ as the drive becomes large ${\lambda }_0\to \mathrm{\infty }$. In a spiking neuron, this corresponds to regular spiking with an inter-spike interval of $\mathrm{\Delta }t=1$. The large-input limit with absolute refractoriness therefore yields a regular spiking limit of the Poisson GLM spiking model (Figure 1c).

Sampling from this model

Sampling from this model with a finite step-size differs slightly from what one would ordinarily do with an autoregressive Poisson process. Ordinarily, one samples $y\sim \mathrm{Poisson}(\lambda \cdot \mathrm{\Delta }t)$. However, because of refractoriness, the probability of observing counts $y>1$ is zero, for $\mathrm{\Delta }t$ shorter than the refractory period.

Instead, we should ask what the probability of not spiking is, within time window $\mathrm{\Delta }t$. The probability of a zero count in a Poisson process is $\exp (-\lambda \cdot \mathrm{\Delta }t)$. The probability of emitting a spike is therefore $p=1-\exp (-\lambda \cdot \mathrm{\Delta }t)$. Equivalently, one may sample Poisson counts, but replace any multiple-spike events with single spikes. Note that this nonlinearity is similar, but not identical, to the sigmoidal nonlinearity that one would use for a Bernoulli GLM

This modified sampling is necessary when treating approximations of the regular spiking limit, since instantaneous conditional intensities diverge toward infinity. The problem is reminiscent of the issue of infinite populations in infinitesimal volume encountered when taking limits in neural field theories, and is solved by explicitly addressing considering activity of a window with a defined volume in time $\mathrm{\Delta }t$.

It's also important to account for how refractoriness changes the variance of this process. Ordinarily, samples of rate $\lambda$ and time-step $\mathrm{\Delta }t$ can be thought of as Poisson with $\mu ={\sigma }^2=\lambda \cdot \mathrm{\Delta }t$. However, the strong inhibition associated with the absolute refractoriness reduces this variance. The variance can instead be modeled as Bernoulli, with ${\sigma }^2=p(1-p)$.

Figure 1: Sampling from a refractory-Poisson autoregressive GLM. (a) For a constant input ${\lambda }_0$, the mean-rate of a refractory Poisson model with an absolute refractory period of $\tau =1$ behaves like a continuous-time version of a Bernoulli process. The output firing rate is a logistic sigmoid function of the log-input $\ln {\lambda }_0$. (b) The sampling and observation model for the refractory-Poisson should take into account both the continuous time nature, and the limited firing rate. The probability of observing a spike is given not by a Bernoulli or Poisson mode, but rather $Pr(y=1)=1-\exp [-\lambda \cdot \mathrm{\Delta }t]$. (c) For low rates, the refractory-Poisson model still exhibits some randomness. Top: spiking output from twenty random samples of the refractory-Poisson model driven by input ${\lambda }_0=10$. Bottom: spiking output for strong drive, ${\lambda }_0=100$, showing more regular spiking.

Continuous approximations

In Rule and Sanguinetti (2018) we used a Langevin approximation to autoregressive point-processes, which approximates the Poisson noise as Gaussian:

$$y=dN=\lambda \cdot dt+\sqrt{\lambda }\cdot dW,$$

where $dN$ denotes the derivative of the count process (the integral of $y(t)$ in time), and $dW$ denotes the derivative of the standard Wiener process. Normally, one integrates the resulting stochastic differential equation using Euler-Maruyama.

$$\begin{array}{rl}\mathrm{\Delta }N& =\lambda \cdot \mathrm{\Delta }t+\sqrt{\lambda \cdot \mathrm{\Delta }t}\cdot \xi \\ \xi & \sim \mathcal{N}(0,1)\end{array}$$

However, this diffusive approximation is wrong in the case of the refractory Poisson process. The refractory period reduces variance, and spiking behaves more like a Bernoulli variable where ${\sigma }^2=p(1-p)$, where $p=1-\exp (-\lambda \cdot \mathrm{\Delta }t)$. For any finite $\mathrm{\Delta }t$, we instead get the integration update:

$$\begin{array}{rl}\mathrm{\Delta }N& =[1-e^{-\lambda \mathrm{\Delta }t}]+\sqrt{p(1-p)}\cdot \xi \\ & =[1-e^{-\lambda \mathrm{\Delta }t}]+e^{-\lambda \mathrm{\Delta }t}\sqrt{e^{\lambda \mathrm{\Delta }t}-1}\cdot \xi \end{array}$$

This reduces to the original approximation when $\lambda \mathrm{\Delta }t\to 0$ is small, but we cannot assume that in the case of the deterministic limit, since the instantaneous rate $\lambda$ can diverge. To get a process modeling regular spiking, we simultaneously take the limits ${\lambda }_0\to \mathrm{\infty }$ and $\mathrm{\Delta }t\to 0$. This does not behave like an Itô stochastic differential equation, but the process is well-defined for all finite $\mathrm{\Delta }t$, and can be numerically integrated.

(or the Complex Chemical Langevin Equation (CCLE) of Schnoerr, Sanguinetti, and Grima, 2014)

Moment equations

In " Moment equations for an autoregressive point process with an absolute refractory period " (February 5, 2019; online notes ), we derived equations for how the mean and covariance evolve in an AR-PPGLM with absolute refractoriness. The moment equations are much simpler for a neuron with steady-state drive, and for which all spiking history effects are captured by an absolute refractory period. For fixed input ${\lambda }_0$ and in the absence of history effects ($\eta (\tau )=0$), the refractory Poisson model is linear. The time evolution of the history process $h(t,\tau )=y(t-\tau )$ can be written as:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\partial }_{t}h& =-{\partial }_{\tau }h+{\delta }_{\tau =0}y\\ y& \sim \mathrm{Poisson}(\lambda \cdot dt)\\ \lambda & =(1-\text{R}h)\cdot {\lambda }_0\end{array}\end{array}$$

Because the interaction between the firing rate and the spiking history is linear, evolution of the moment (mean $\mu$ and covariance $\mathrm{\Sigma }$) of the process history can be written without invoking any approximations or moment closures needed:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mu }_y& =(1-\text{R}{\mu }_h){\lambda }_0\cdot dt\\ {\partial }_t\mu & =-{\partial }_{\tau }\mu +{\delta }_{\tau =0}{\mu }_y\\ {\partial }_t\mathrm{\Sigma }& =J\mathrm{\Sigma }+\mathrm{\Sigma }{J}^{\mathrm{\top }}+{\delta }_{\tau =0}{\mu }_y{\delta }_{\tau =0}^{\mathrm{\top }}\\ J& =-{\partial }_{\tau }-{\delta }_{\tau =0}{\text{R}}^{\mathrm{\top }}{\lambda }_0\cdot dt\end{array}\end{array}$$

Figure 2 illustrates that these moment equations match the true process, and also checks the accuracy of the continuous approximation derived earlier.

Figure 2: Solution for the mean and covariance of the refractory Poisson model. In this simulation, the drive is ${\lambda }_0=25$, the absolute refractory period is $\tau =1$, and time step is $\mathrm{\Delta }t=0.02$. (a) Single-time marginal means and variances. Black and shaded region trace represents the theoretical mean and variance, computed by integrating the moment equation. Green traces reflect $\pm 1\sigma$ computed by sampling the refractory Poisson model. Teal traces reflect $\pm 1\sigma$ sampled from a continuous approximation, with a Bernoulli-like, rather than Poisson, noise model. (b) Similar to (a), but in this case means and variances have been averaged over one refractory period. The refractory period reduces the variability of the process, so the average fluctuations over a short time period are smaller than one might expect given the single-time marginal variances of the process.

Overall,

These notes examined one way to model an absolute refractory period, and its role in reducing neuronal variability, in GLM point-process models. We showed that the sampling of the process must be modified to account for the influence of the refractory period, which enforces a sort of "law of conservation of spikes" in time. The Langevin approximation to the refractory Poisson model is also not straightforward, but for any finite $\mathrm{\Delta }t$ one can derive an integration scheme that approximates the original process. The resulting process resembled a continuous-time Bernoulli process, but with slightly different nonlinearities. In the case of steady-state input, the moment equations for the refractory Poisson model are linear, and agree with Monte-Carlo simulations. It might be interesting to see whether other types of highly regular spiking observed in neurons, such as integrators or resonators with variable frequencies, can be defined in a consistent way in terms of AR-PPGLMs.