Moment equations for an autoregressive point process with an absolute refractory period

These notes extend the autoregressive point-process moment closures derived in Rule and Sanguinetti (2018) to models that include an absolute refractory period via a gating term that sets the rate to zero after a spike.

[get notes as PDF]

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,\tau )+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(t,\tau )=\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,\tau )]\cdot \varphi (\text{H}h(t,\tau )+u(t))\end{array}\end{array}$$

The gating function $g[h(t,\tau )]$ is defined as:

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

The coefficients $\rho (\tau )$ are a refractory history filter. For example, to constrain a neuron to fire at most once per unit time, set $\rho (\tau )=1$ for $\tau \in (0,1]$ (and zero otherwise). If we define a linear operator $R$ as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\text{R}h(t,\tau )& =\underset{ϵ\downarrow 0}{lim}{\int }_{ϵ}^1\rho (\tau )h(t,\tau )d\tau ,\end{array}\end{array}$$

then we may write the firing rate as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\lambda & =(1-\text{R}h)\cdot \varphi (\text{H}h+u).\end{array}\end{array}$$

In this notation, the history process $h$ evolves as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\partial }_{t}h& =-{\partial }_{\tau }h+{\delta }_{\tau =0}y\\ y& \sim \mathrm{Poisson}(\lambda \cdot \mathrm{\Delta }t)\\ \lambda & =(1-\text{R}h)\cdot \varphi (\text{H}h+u)\end{array}\end{array}$$

We are interested in deriving equations for the mean and covariance of the auxiliary history process $h$. To incorporate corrections from covariance effects, we need to estimate the expected (mean) rate $⟨\lambda ⟩$, based on the mean and covariance of the spiking history:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}⟨\lambda ⟩& =⟨(1-Rh)\varphi (Hh+u)⟩.\end{array}\end{array}$$

We'll the a Gaussian moment closure on a local, second-order Taylor expansion of the (nonlinear) intensity described in Rule and Sanguinetti (2018).

If we assume that $h$ is Gaussian, then correcting for higher moments amounts to correcting for the interaction between the covariance $\mathrm{\Sigma }$ and the curvature of $\lambda$.

This moment closure is derived for models without absolute refractoriness, $\lambda =\varphi (\text{H}h+u)$, in Rule and Sanguinetti (2018). In these notes, we therefore focus on the additional term $\text{R}h\cdot \varphi (\text{H}h+u)$ arising from the absolute refractory period.

Covariance corrections to the mean

To derive a covariance correction to the mean, we need to calculate the curvature of the term $\text{R}h\varphi (\text{H}h+u)$ about the current mean $\mu (t,\tau )=⟨h(t,\tau )⟩$.

In the following derivations, we denote functions of time-lag $\tau$ $h({\tau }_1,t)$, $h({\tau }_2,t)$ as $h_i$, $h_j$, etc., for compactness. Using the abbreviations

$$\begin{array}{rl}\varphi & =\varphi (\text{H}h+u),\\ {\varphi }^{\prime }& ={\varphi }^{\prime }(\text{H}h+u),\end{array}$$

etc. The derivatives are

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\partial }_{h_i}[\text{R}h\varphi ]& =(\text{R}h){\varphi }^{\prime }{\eta }_i+{\rho }_i\varphi \\ {\partial }_{h_j}{\partial }_{h_i}[\text{R}h\varphi ]& ={\rho }_j{\varphi }^{\prime }{\eta }_i+(\text{R}h){\varphi }^{″}{\eta }_j{\eta }_i+{\rho }_i{\varphi }^{\prime }{\eta }_j\end{array}\end{array}$$

Or in vector notation:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{\nabla }\left[\text{R}h\varphi \right]& =(\text{R}h){\varphi }^{\prime }\circ \text{H}+\varphi \circ \text{R}\\ \mathrm{\nabla }{\mathrm{\nabla }}^{\mathrm{\top }}\left[\text{R}h\varphi \right]& ={\varphi }^{\prime }\circ [\text{R}{\text{H}}^{\mathrm{\top }}+\text{H}{\text{R}}^{\mathrm{\top }}]+(\text{R}h){\varphi }^{″}\circ \text{H}{\text{H}}^{\mathrm{\top }},\end{array}\end{array}$$

where $\circ$ denotes element-wise multiplication.

The covariance corrections to the mean come from how the noise (covariance) interacts with the second derivatives derived above:

$$\frac{1}{2}\sum _{ij}\mathrm{\Delta }{h}_i\mathrm{\Delta }{h}_j\cdot {\partial }_{h_j}{\partial }_{h_i}[\dots ].$$

In vector notation, this can be abbreviated as

$$\frac{1}{2}\mathrm{tr}[\mathrm{\Sigma }\cdot \mathrm{\nabla }{\mathrm{\nabla }}^{\mathrm{\top }}(\dots )].$$

The covariance correction without absolute refractoriness is (Rule and Sanguinetti, 2018):

$$\begin{array}{}\text{(10)}& \begin{array}{r}\frac{1}{2}{\varphi }^{″}\circ \text{H}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }},\end{array}\end{array}$$

and the covariance corrections for the absolute-refractory term $⟨(\text{R}h)\varphi ⟩$ are

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \frac{1}{2}\mathrm{tr}[\mathrm{\Sigma }\cdot ({\varphi }^{\prime }\circ [\text{R}{\text{H}}^{\mathrm{\top }}+\text{H}{\text{R}}^{\mathrm{\top }}]+(\text{R}\mu ){\varphi }^{″}\circ \text{H}{\text{H}}^{\mathrm{\top }})]\\ =& \frac{1}{2}\{{\varphi }^{\prime }\circ [\text{R}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}+\text{H}\mathrm{\Sigma }{\text{R}}^{\mathrm{\top }}]+(\text{R}\mu ){\varphi }^{″}\circ \text{H}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}\}.\end{array}\end{array}$$

Combining these, we get the overall correction to the mean of

$$\begin{array}{}\text{(12)}& \begin{array}{r}\frac{1}{2}[(1-\text{R}\mu ){\varphi }^{″}\circ \text{H}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}+{\varphi }^{\prime }\circ (\text{R}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}+\text{H}\mathrm{\Sigma }{\text{R}}^{\mathrm{\top }})].\end{array}\end{array}$$

I personally find the notation a bit more readable if we define the variance of the input ${\sigma }_a^2=\text{H}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}$, the gating term $g=1-\text{R}\mu$, and a variable $K=\text{R}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}$ that summarizes interactions between the refractory gating and autohistory filter:

$$\begin{array}{}\text{(13)}& \begin{array}{r}\frac{1}{2}[g{\varphi }^{″}\circ {\sigma }_a^2+{\varphi }^{\prime }\circ (K+K^{\mathrm{\top }})]\end{array}\end{array}$$

In words: the variance in the activation, ${\sigma }_a^2=\text{H}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}$, interacts with the curvature of the firing rate nonlinearity ${\varphi }^{″}$. This is similar to the model without the absolute refractory period, just multiplied by the mean gating term $g=1-\text{R}\mu$. We also get an additional correction from the slope of the nonlinearity, ${\varphi }^{\prime }$. This is related to how terms in the absolute refractory filter $\text{R}$ correlate with terms in the history filter $\text{H}$, $K=\text{R}\mathrm{\Sigma }{\text{H}}^{\mathrm{\top }}$.

Time evolution of the covariance

Following Rule and Sanguinetti (2018), the covariance evolves according to the the Jacobian of the mean system, similar to the extended-time Kalman-Bucy filter. Refer to Rule and Sanguinetti (2018) for the linear terms. The introduction of an absolute refractory filter only changes the nonlinear terms in the Jacobian, $\mathrm{\nabla }\lambda$. This is given by:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathrm{\nabla }\lambda & =g{\varphi }^{\prime }\circ \text{H}-\varphi \circ \text{R}.\end{array}\end{array}$$

These moment equations generalize to population models, where the history and refractory filters are replaced by matrices of filters for each member of the population, and history moments now refer to the joint moments of the population history.