Multitaper Spectral Analysis of Neuronal Spiking Activity Driven by Latent Stationary Processes

TL;DR

This paper introduces a specialized multitaper spectral estimation method for neuronal spiking data, enabling more accurate analysis of brain rhythms underlying neural activity.

Contribution

It develops a novel variant of the multitaper method tailored for point process data, improving spectral estimation accuracy in neuroscience applications.

Findings

Significant reduction in bias and variance compared to existing methods

Efficient computation of spectral estimates from spiking data

Enhanced ability to infer latent neural processes

Abstract

Investigating the spectral properties of the neural covariates that underlie spiking activity is an important problem in systems neuroscience, as it allows to study the role of brain rhythms in cognitive functions. While the spectral estimation of continuous time-series is a well-established domain, computing the spectral representation of these neural covariates from spiking data sets forth various challenges due to the intrinsic non-linearities involved. In this paper, we address this problem by proposing a variant of the multitaper method specifically tailored for point process data. To this end, we construct auxiliary spiking statistics from which the eigen-spectra of the underlying latent process can be directly inferred using maximum likelihood estimation, and thereby the multitaper estimate can be efficiently computed. Comparison of our proposed technique to existing methods using simulated data reveals significant gains in terms of the bias-variance trade-off.

Tables

Table 1. Table 1: Normalized MSE Comparisons

PMTM	SS-PSD	PSTH-PSD
$0.4733\pm 0.0072$	$0.8164\pm 7.9592\times {10}^{-8}$	$7.7772\pm 2.0641$

Equations

$$\lambda (t|H_t):=\underset{\Delta \to 0}{\lim }\frac{P[N(t+\Delta )-N(t)=1|H_t]}{\Delta }.$$

$$x_k={\int }_{-\frac{1}{2}}^{\frac{1}{2}}{e}^{i2\pi fk}dz(f),$$

$${\lambda }_k=\mu +x_k,n_k^{(l)}\sim \text{Bernoulli}({\lambda }_k),$$

$$\displaystyle x_{k}\!=\!\sum_{m=1}^{N}\frac{2}{N}\Big{[}a_{m}\cos\frac{\pi(m\!-\!1)}{N}-b_{m}\sin\frac{\pi(m\!-\!1)}{N}\Big{]},$$

$$\mathbf{A}:=\frac{2}{N}\left[\begin{array}{cccccc}1& \cos \frac{\pi }{N}& -\sin \frac{\pi }{N}& \dots & \cos \frac{(N-1)\pi }{N}& -\sin \frac{(N-1)\pi }{N}\\ 1& \cos \frac{2\pi }{N}& -\sin \frac{2\pi }{N}& \dots & \cos \frac{2(N-1)\pi }{N}& -\sin \frac{2(N-1)\pi }{N}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& \cos \frac{K\pi }{N}& -\sin \frac{K\pi }{N}& \dots & \cos \frac{K(N-1)\pi }{N}& -\sin \frac{K(N-1)\pi }{N}\end{array}\right],$$

$$\frac{f_i(z_i)\text{mathbbm}1[-\mu \le z_i\le \mu ]}{\int f_i(\xi )\text{mathbbm}1[-\mu \le \xi \le \mu ]d\xi },$$

$${\hat{S}}^{(j)}(f):=$$

$${\hat{S}}^{(\mathsf{mtm})}(f):=\frac{1}{J}\sum _{j=1}^J{\hat{S}}^{(j)}(f).$$

$$n_k^{(l,j)}=n_k^{(l)}{v}_k^{(j)}\text{mathbbm}1[v_k^{(j)}\ge 0]-{\check{n}}_k^{(l)}{v}_k^{(j)}\text{mathbbm}1[v_k^{(j)}<0]$$

$${\lambda }_k^{(j)}={\mu }_k^{(j)}+x_k{v}_k^{(j)},$$

$${\hat{\mathbf{\theta }}}_{\mathsf{ML}}^{(j)}=\underset{{\mathbf{\theta }}^{(j)}}{\mathrm{argmax}}P({\mathcal{D}}^{(j)}|{\mathbf{\theta }}^{(j)})$$

$$\log L(\mathbf{\theta }|\mathbf{z},\mathcal{D})=\sum _{l=1}^L$$

$$\displaystyle-\sum_{m=1}^{2N-1}\Big{(}\log\int_{\!-\!\mu}^{\mu}\!f_{m}(\mathbf{\xi})d\mathbf{\xi}+\frac{z_{m}^{2}}{2\sigma_{m}^{2}}+\frac{1}{2}\log\sigma_{m}^{2}\Big{)}+{C},$$

$$\displaystyle\mathbf{Q}(\boldsymbol{\theta}|\boldsymbol{\theta}^{[i]})=-\sum_{m=1}^{2N-1}\bigg{(}\frac{1}{2}\log\sigma_{m}^{2}+\frac{1}{2\sigma_{m}^{2}}\mathbbm{E}[z_{m}^{2}|\mathcal{D},\boldsymbol{\theta}^{[i]}]\bigg{)}+{C^{\prime}},$$

$$\displaystyle\boldsymbol{\mu}_{\mathbf{z}^{[i]}}=\operatorname*{arg\,max}_{\mathbf{z}\in D}\ \ \sum_{l=1}^{L}\sum_{k=1}^{K}\bigg{[}n_{k}^{(l)}\log\frac{\mu_{k}+(\mathbf{Az})_{k}}{1-(\mu_{k}+(\mathbf{Az})_{k})}+\log{(1\!-\!(\mu_{k}+(\mathbf{Az})_{k}))}\bigg{]}-\sum_{m=1}^{2N-1}\frac{z_{m}^{2}}{2\sigma_{m}^{2}}\ ,$$

Full text

1 Introduction

Spectral analysis techniques are among the most important tools for extracting information from time-series data recorded from naturally occurring processes. Examples include speech [1], image and video [2], and electroencephalography (EEG) [3] data. Due to the exploratory nature of most of these applications, nonparametric techniques based on Fourier methods and Wavelets are among the most widely used. In particular, the multitaper method excels among the available nonparametric techniques due to both its simplicity and control over the bias-variance trade-off by means of bandwidth adjustment[4, 5, 6]. This technique has been successfully utilized in the analyze of EEG data [7, 8, 6].

The advent of large-scale invasive recording technologies from the brains of animals and humans, such as multi-electrode arrays and electrocorticography (ECoG), has resulted in abundant pools of neuronal spiking data, which often exhibit oscillatory features [9]. Characterizing the properties of these oscillations with high spectral resolution is crucial to understanding their role in cognitive functions.

Most existing spectral analysis techniques, however, are designed for continuous-time data and cannot be readily applied to binary spiking data (See, for example, [10] and [11]). There have been efforts aimed at addressing this challenge, which consider the periodogram of smoothed spike trains using kernel methods as the spectral representation [12, 13, 14]. Another strand of results are based on the theory of point processes, which has been widely used in recent years to model and analyze the statistical properties of binary spike trains [15, 16, 17, 18]. These techniques relate the Conditional Intensity Function (CIF) or the spiking rate of a point process governing the spiking statistics to intrinsic and external neural covariates using state-space models. Then, the spectrum of the estimated CIF is characterized using standard nonparametric or parametric techniques[19, 20, 21].

Despite their relative success in application, these methods have several shortcomings from a theoretical perspective. First, it is known that direct smoothing of the signal using kernels or indirect smoothing using state-space models, results in the distortion of the spectrum [22]. Second, time-domain smoothing alleviates the variance of the estimates at the cost of increasing the bias. On the other hand, existing techniques which avoid time-domain smoothing (e.g., [23]) may exhibit high variability. Third, these modeling frameworks often require a priori information (e.g., sparsity) or may suffer from model mismatch (e.g., overly-smoothed state estimates).

To address these issues, in this paper we introduce a novel multitaper spectral analysis method, which we call the Point Process Multitaper Method (PMTM), to be directly applied to binary data. To this end, we generate auxiliary spiking statistics which correspond to the tapered versions of the CIF, which are then used to independently estimate the eigen-spectra of the tapered CIFs via the Maximum Likelihood (ML) procedure. The multitaper spectral estimate is formed by averaging the corresponding eigen-spectral estimates. Our approach distinguishes itself from existing work by directly estimating the spectra from binary observations via a novel adaptation of multitaper analysis, with no recourse to smoothing or need for a priori information. We demonstrate the performance of PMTM using simulated spike trains driven by an autoregressive (AR) process. Our results reveal substantial gains achieved by PMTM as compared to existing nonparametric techniques, in terms of the bias-variance trade-off.

2 Problem Formulation

Let $N(t)$ and $H_{t}$ denote the point process representing the number of spikes and spiking history in $[0,t)$, respectively, where $t\in[0,T]$ and $T$ denotes the observation duration. The CIF of a point process $N(t)$ is defined as:

[TABLE]

To discretize the continuous process, we consider time bins of length $\Delta$, small enough that the probability of having two or more spikes in an interval of length $\Delta$ is negligible. Thus, the discretized point process can be modeled by a Bernoulli process with success probability $\lambda_{k}:=\lambda(k\Delta|H_{k})\Delta$, for $1\leq k\leq K$, where $K:=T/\Delta$ and is assumed to be an integer with no loss of generality. Note that $\lambda_{k}$ forms the CIF of the discretized process, which we refer to as CIF hereafter for brevity. Let $n_{k}\in\{0,1\}$ be the number of spikes in bin $k$, for $0\leq k\leq K$. Our objective is to estimate the Power Spectral Density (PSD) of the CIF from the observed spike train $\{n_{k}\}_{k=1}^{K}$, under the assumption that the CIF is a second-order stationary process.

More generally, we consider an ensemble of $L$ neurons or $L$ trials from a single neuron driven by the same CIF, and denote the observed spike trains by \mathcal{D}:=\big{\{}n_{k}^{(l)}\big{\}}_{k=1,l=1}^{K,L}. When considering an ensemble of $L$ neurons, this setting may only be valid for neuronal recordings from a small area of cortex using multi-electrode arrays (See, for example [24]), and when considering $L$ trials from the same neuron, it is assumed that all trials pertain to the same stimulus (See, for example [25]). Nevertheless, in many other cases of interest, one has only access to a single trial from a single neuron, which makes spectral estimation an even more challenging problem. We model the CIF using a second order zero-mean stationary random process, $\{x_{k}\}_{k=1}^{K}$, which by virtue of Spectral Representation Theorem [22] admits a Cramér representation [26] of the form:

[TABLE]

where $dz(f)$ is a complex-valued orthogonal increment process and the PSD, $S(f)$ of the process is defined as: $S(f)df=\mathbbm{E}[|dz(f)|^{2}]$. Finally, we use a linear link for the CIF so that model can be summarized as

[TABLE]

for $\ 1\leq k\leq K$ and $1\leq l\leq L$. The choice of the linear link, as opposed to more natural links such as the logistic function, is for the sake of simplicity of the auxiliary data generation process that will be described in Section 3.1. Acknowledging the non-linearity of the model (3) and availability of only a finite number of samples, we consider a piece-wise continuous approximation to the PSD, i.e, $dz(f)$ is constant over the intervals $[\frac{m-1}{2N},\frac{m}{2N})$, for large enough $N$, for $m=1,2,\cdots,N$. This enables us to express $dz(f)=(a_{m}+ib_{m})df$ for $f\in[\frac{m-1}{2N},\frac{m}{2N})$, where $a_{m}$ and $b_{m}$ are random variables for $m=1,2,\cdots,N$ [23]. Invoking the conjugate symmetry of $dz(f)$ for real valued $\{x_{k}\}_{k=1}^{K}$, Eq. (2) can be written as

[TABLE]

with a PSD of $S(f)=\frac{1}{N}\mathbbm{E}[a_{m}^{2}+b_{m}^{2}]$ for $f\in[\frac{m-1}{2N},\frac{m}{2N})$.

Denoting $\mathbf{x}:=[x_{1},x_{2},\cdots,x_{K}]^{\top}$ and $\mathbf{z}:=[a_{1},a_{2},b_{2},\cdots,a_{N},b_{N}]^{\top}$ and defining $\mathbf{A}$ as

[TABLE]

one can write (4) in the vector form as $\mathbf{x}=\mathbf{Az}$. By the orthogonality of the increment process $dz(f)$, $z_{i}$’s are independent. We further assume that $z_{i}$, for $i=1,2,\cdots,2N-1$, follows a truncated zero-mean Gaussian distribution, with a density

[TABLE]

where $f_{i}(\cdot)$ is the Gaussian density $\mathcal{N}(0,{\sigma_{i}}^{2})$ and $\mathbbm{1}$ is the indicator function. This ensures that each entry of $\{x_{k}\}_{k=1}^{K}$ is restricted to $[-\mu,\mu]$ for some positive scalar $\mu$ and by selecting $\mu\leq 1/2$, one can achieve $0\leq\lambda_{k}\leq 1$, for all $1\leq k\leq K$.

3 The Point Process Multitaper Method

The MTM is an extension of tapered PSD estimation, where the spectral estimate is computed by averaging several tapered PSD estimates corresponding to orthogonal tapers with optimal spectral leakage properties [22]. The set of tapers from the Discrete Prolate Spheroidal Sequences (dpss) [27] provides excellent control over the bias-variance trade-off [6, 4, 5].

Let $v_{k}^{j}$ is be the $k^{\text{th}}$ sample of the $j^{\text{th}}$ dpss sequence, for a given design bandwidth $W$, $k=1,2,\cdots,K$ and $j=1,2,\cdots,J$ such that $J<\lfloor 2KW\rfloor-1$. Given the time-series data \big{\{}x_{k}\big{\}}_{k=1}^{K}, the $j^{\text{th}}$ eigen-spectrum is given by:

[TABLE]

from which the MTM PSD estimate can be computed as:

[TABLE]

Due to the non-linear nature of the model (3), forming the MTM PSD estimate based on spiking data becomes non-trivial. In what follows, we indeed address this issue by devising a novel variant of MTM.

3.1 Generating Auxiliary Spiking Statistics

Given that $\{\lambda_{k}\}_{k=1}^{K}$ is not directly observable, we instead generate auxiliary spiking statistics which would have been generated by the tapered CIFs, i.e., \big{\{}v_{k}^{(j)}\lambda_{k}\big{\}}_{k=1}^{K}. For non-negative $v_{k}^{(j)}$ (e.g., for $j=1$), generating such modified spike trains is usually carried out using the thinning method [28]. The basic idea of the thinning method is to retain the original spikes with a probability determined by the ratio of the target and original CIFs. Noting that the ensemble average of the spike trains form a sequence of sufficient statistics for estimating the CIF, one can generate several realizations of thinned spike trains and consider their ensemble average as a sufficient statistic for the thinned process. It is not difficult to see that for non-negative tapers, the ensemble average converges in probability to the original spiking activity multiplied element-wise by the taper [29]. Thus, for the sake of robustness we use this limit as the ensemble average of the auxiliary spike trains.

Given that the dpss tapers take negative values (for $j>1$), the aforementioned thinning method for constructing spike trains from a tapered CIF cannot be applied, as its naïve application results in negative-valued spikes. To resolve this issue, we use the fact that the spike train is generated according to a Bernoulli process and therefore its complement given by $\check{n}^{(l)}_{k}:=1-n^{(l)}_{k}$ has a CIF given by $\check{\lambda}_{k}:=1-\lambda_{k}=1-\mu-x_{k}$. Therefore, the non-negative sequence

[TABLE]

has a limit ensemble average corresponding to the CIF:

[TABLE]

where $\mu_{k}^{(j)}=\mu v^{(j)}_{k}\mathbbm{1}[v^{(j)}_{k}\!\geq\!0]\!-\!(1\!-\!\mu)v^{(j)}_{k}\mathbbm{1}[v^{(j)}_{k}\!<\!0]$, for $j=1,2,\cdots,J$, and can be utilized as a sequence of sufficient statistics to estimate the spectral representation of the tapered process. Note that the non-negativity of the sequence $n_{k}^{(l,j)}$ follows from the definition of $\check{n}_{k}^{(l)}$. In addition, each taper is scaled by its maximum absolute value to ensure $n_{k}^{(l,j)}\leq 1$, and the estimated eigen-spectra are accordingly rescaled. Using this approach, we can generate the ensemble average of auxiliary spike trains for any taper regardless of its sign. It is noteworthy that while in principle this procedure can be extended to more general link functions, the generation of the corresponding auxiliary statistics may be more intricate. Fig. 1 provides a visual summary of our proposed framework. The time-bandwidth product and the number of tapers are chosen following guidelines from the MTM literature [4, 6].

3.2 Maximum Likelihood Estimation of the Eigen-spectra

Once the auxiliary spiking statistics \mathcal{D}^{(j)}=\big{\{}n^{(l,j)}_{k}\big{\}}_{k=1,l=1}^{K,L}, $j=1,2,\cdots,J$ are available, the eigen-spectra need to be estimated to construct the PSD. Given the modeling framework of Section 2, estimation of the $j^{\text{th}}$ eigen-spectrum reduces to estimating the parameters $\boldsymbol{\theta}^{(j)}:=[{\sigma_{1}^{(j)2}},{\sigma_{2}^{(j)2}},\cdots,{\sigma^{(j)2}_{2N-1}}]^{\top}$, where ${\sigma_{m}^{(j)2}}$ is the variance of the random variable $z_{i}^{(j)}$ corresponding to the $j^{\text{th}}$ eigen-spectra, $i=1,2,\cdots,2N-1$. The ML estimate of the parameter $\boldsymbol{\theta}^{(j)}$ is given by:

[TABLE]

Note that expressing $P(\mathcal{D}^{(j)}|\boldsymbol{\theta}^{(j)})$ solely in terms of $\mathcal{D}^{(j)}$, i.e., eliminating $\mathbf{z}^{(j)}:=[z_{1}^{(j)},z_{2}^{(j)},\cdots,z_{2N-1}^{(j)}]^{\top}$, introduces computational intricacies, which we avoid by using the Expectation-Maximization (EM) algorithm [30] as our solution method. In what follows we drop the superscript $j$ for the sake of clarity, as the same procedure will be used for estimating each eigen-spectrum. If $\mathbf{z}$ is known, the complete data log-likelihood of the observations $\mathcal{D}$ can be written as:

[TABLE]

which could be efficiently maximized to estimate $\boldsymbol{\theta}$ (terms independent of $\boldsymbol{\theta}$ are denoted by $C$).

Before presenting the EM algorithm, we note that $\int_{-\mu}^{\mu}f_{m}({\xi})d{\xi}=1-2{\Phi}(\gamma_{m})\approx 1$ for moderate values of $\gamma_{m}:=\frac{\mu}{\sigma_{m}}$ and small enough $\sigma_{m}^{2}$. Thus, we drop this term henceforth to avoid unnecessary complexity. One may choose to work with this term included, at the expense of additional computational costs. At the $i^{\text{th}}$ iteration, we have:

3.2.1 E-step

Given $\boldsymbol{\theta}^{[i]}$, the Q-function is given by

[TABLE]

which requires $f_{\mathbf{z}|\mathcal{D},\boldsymbol{\theta}^{[i]}}$ or samples from it, and can thus be computationally demanding to compute (terms independent of $\boldsymbol{\theta}$ are denoted by $C^{\prime}$). Instead, we use the unimodality of the density and approximate it by a multivariate Gaussian density $\mathcal{N}(\boldsymbol{\mu}_{\mathbf{z}^{[i]}},\boldsymbol{\Sigma}_{\mathbf{z}^{[i]}})$ [19, 15]. By invoking the fact that the mode and mean of a multivariate Gaussian density coincide and the Hessian of its natural logarithm is equal to -\big{(}{\boldsymbol{\Sigma}_{\mathbf{z}^{[i]}}\big{)}}^{-1}, we get:

[TABLE]

and $\boldsymbol{\Sigma}_{\mathbf{z}^{[i]}}$ is given by the Hessian of the log-likelihood in (3.2) evaluated at $\boldsymbol{\mu}_{\mathbf{z}^{[i]}}$. The maximization problem (13) is concave over $D=\{\mathbf{z}\in\mathbb{R}^{2N-1}:{0}\leq\mu_{k}+(\mathbf{Az})_{k}\leq{1},k=1,2,\cdots,K\}$ and the Hessian is negative definite, so Newton-type method for bound-constrained optimization can be used to compute $\boldsymbol{\mu}_{\mathbf{z}^{[i]}}$ efficiently. We use a line-search method [31], which generates a sequence of iterates by setting $\boldsymbol{\mu}^{[r+1]}_{\mathbf{z}^{[i]}}=\boldsymbol{\mu}^{[r]}_{\mathbf{z}^{[i]}}+\alpha^{[r]}\mathbf{d}^{[r]}$, where $\boldsymbol{\mu}^{[r+1]}_{\mathbf{z}^{[i]}}$ is a feasible approximation to the solution, $\alpha^{[r]}$ is the step-size and $\mathbf{d}^{[r]}$ is the Newton’s step for that iteration. Then, $\boldsymbol{\Sigma}_{\mathbf{z}^{[i]}}$ can be computed by evaluating the Hessian of (3.2) at $\mathbf{z}=\boldsymbol{\mu}_{\mathbf{z}^{[i]}}$, which allows $\mathbbm{E}[z_{m}^{2}|\mathcal{D},\boldsymbol{\theta}^{[i]}]$ to be calculated as \big{(}(\boldsymbol{\mu}_{\mathbf{z}^{[i]}})_{m}\big{)}^{2}+\big{(}\boldsymbol{\Sigma}_{\mathbf{z}}^{[i]}\big{)}_{m,m}.

3.2.2 M-step

The parameter vector $\boldsymbol{\theta}^{[i+1]}$ is updated by maximizing the expectation in (12). Given that $Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{[i]})$ is concave over the positive orthant, its unique maximizer is given by $\widehat{\boldsymbol{\theta}}_{m}^{[i+1]}=\mathbbm{E}[z_{m}^{2}|\mathcal{D},\boldsymbol{\theta}^{[i]}]$.

Note that we have assumed $\mu_{k}^{(j)}$’s to be known. Since it is not the case for most practical purposes, we first estimate $\mu$ as $\hat{\mu}=\frac{1}{LK}\sum_{l,k=1}^{L,K}n_{k}^{(l)}$ and compute $\mu_{k}^{(j)}$ in (9) using $\hat{\mu}$. We terminate the EM algorithm after a fixed large number of iterations or until some convergence criterion is met. A similar stopping rule for the maximization problem inside each EM step is used. We initialize $\boldsymbol{\theta}^{[0]}$ as an arbitrary vector in the positive orthant. Following the termination of the EM algorithm, the eigen-spectra are calculated as $\widehat{S}(0)=\widehat{\sigma}_{1}^{2}$ and $\widehat{S}(f_{m})=\widehat{\sigma}_{2m}^{2}+\widehat{\sigma}_{2m+1}^{2}$ for $f_{m}=\frac{m}{2N}$ and $m=1,2,\cdots N-1$. Finally, the PMTM estimate can be computed using (7). Algorithm 1 summarizes the proposed PMTM procedure.

4 Simulation Results

We simulate $x_{k}$ as an AR$(4)$ process given by $\displaystyle x_{k}=0.4152x_{k-1}-0.0922x_{k-2}+0.4170x_{k-3}-0.8852x_{k-4}+0.025\epsilon_{k}$, where $\epsilon_{j}$ is zero-mean i.i.d. Gaussian noise with unit variance. We compute the CIF as $\lambda_{k}=\mu+x_{k}$ for $\mu=0.12$ (truncated to $[0,1]$, if necessary), to generate the binary spiking activity for $K=512$ samples. A snapshot of one realization of this AR process and the raster plot of $L=10$ spike trains are depicted in Fig. 2A and B, respectively.

We apply PMTM to this simulated data and benchmark it against two existing methods: $(1)$ PSTH-PSD, where the PSD is computed by forming the MTM estimate of the ensemble peristimulus time histogram (PSTH), i.e., the average spike trains, and $(2)$ SS-PSD, where $x_{k}$ is first estimated using a state-space model $x_{k}=x_{k-1}+w_{k}$, followed by forming its MTM PSD estimate [19].

Fig. 3A shows the PMTM (black), PSTH-PSD (green), SS-PSD (aqua) and the true PSD (blue) for the realization shown in 2 in log-scale. For comparison purposes, we have also included the MTM PSD estimate of $x_{k}$, assuming that an oracle has access to it (Oracle PSD, in red). We have used the first $8$ dpss tapers corresponding to $\alpha=5$. As it can be observed from Fig. 3A, PSTH-PSD suffers from high bias, though enjoying reduced variability, and the spectral peaks are difficult to distinguish from the background. On the other hand, SS-PSD suffers from model mismatch, as it over-smooths the CIF due to the usage of a state-space model, and as a result the spectral peaks are nearly absent in the estimate. The PMTM estimate, however, closely follows the true PSD by guarding against spectral leakage and producing a nearly unbiased estimate on par with the Oracle PSD estimate, though it exhibits some variability. In order to quantify these comparisons, we computed a normalized measure of MSE by averaging the squared-error of the PSD normalized by the true PSD values in the log-scale. The normalized MSE values ($\pm$ 2 STD) corresponding to $10$ different AR process realizations and $5$ different spike-train ensemble realizations are given in Table 1, which corroborates our foregoing qualitative comparison. To ease reproducibility, we have deposited a MATLAB implementation of PMTM on the open source repository Github [32], which fully regenerates Fig. 3A.

Finally, Fig. 3B examines the improvement of the PMTM estimates with respect to the ensemble size, for $L=5,10,15$ and $20$. As $L$ increases, the PSD estimates improve, but with a seemingly saturating effect for $L\geq 10$. Note that while the average spiking rate of $\mu=0.12$ is reasonable for some neurons (e.g., Parvalbumin-positive interneurons [33]), often times neurons exhibit spiking rates as low as $\mu=0.05$, for which the performance of all existing methods significantly degrades. Extension of our method for obtaining robust spectral estimates in the domain of low spiking rate and small number of trials is a future direction of research.

5 Conclusion

Spectral estimation of continuous time-series is a well-established domain, as hallmarked by the multitaper method known for its favorable control over the bias-variance trade-off. Computing the spectral representation of the neural covariates that underlie spiking activity, however, sets forth various challenges due to the intrinsic non-linearities involved. In this paper, we addressed this problem by proposing a multitaper method specifically tailored for binary spiking data, which we refer to as PMTM. We compared the performance of PMTM to that of two existing techniques using simulated data, which revealed significant gains in terms of estimation accuracy. The PMTM can be extended to a wide variety of binary data, such as rainfall and earthquake data, to extract spectral representations of the underlying latent processes.

Acknowledgment

This work is supported in part by the National Science Foundation Awards No. 1552946 and 1807216.