Multitaper Analysis of Evolutionary Spectra from Multivariate Spiking Observations

TL;DR

This paper introduces a multitaper spectral estimation method for analyzing the evolutionary spectra of latent non-stationary processes from multivariate neural spiking data, improving bias-variance trade-offs.

Contribution

We develop a novel multitaper spectral estimation technique tailored for multivariate spiking observations, with theoretical bias-variance bounds and demonstrated superior performance.

Findings

Significant bias-variance trade-off improvements over existing methods

Theoretical bounds on estimator bias and variance

Successful application to real neural spiking data

Abstract

Extracting the spectral representations of the neural processes that underlie spiking activity is key to understanding how the brain rhythms mediate cognitive functions. While spectral estimation of continuous time-series is well studied, inferring the spectral representation of latent non-stationary processes based on spiking observations is a challenging problem. In this paper, we address this issue by developing a multitaper spectral estimation methodology that can be directly applied to multivariate spiking observations in order to extract the evolutionary spectral density of the latent non-stationary processes that drive spiking activity, based on point process theory. We establish theoretical bounds on the bias-variance trade-off of the proposed estimator. Finally, we compare the performance of our proposed technique with existing methods using simulation studies and application to real data, which reveal significant gains in terms of the bias-variance trade-off.

Tables

Table 1. TABLE I: Comparison of Relative MSE Performance

Estimation method	Average MSE	Variance of MSE
Oracle ESD	0.0485	$8.2612\times {10}^{-6}$
PPMT-ESD	0.1864	$7.4279\times {10}^{-5}$
SS-ESD	0.3911	$1.3991\times {10}^{-5}$
PSTH-ESD	1.4794	$1.1\times {10}^{-3}$

Equations

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

$$n_{k,j}^{(l)}\sim \mathrm{Bernoulli}({\lambda }_{k,j}),$$

$$x_{k,j}-{\mu }_{k,j}={\int }_{-\pi }^{\pi }{e}^{ik\omega }{A}_{k,j}(\omega )d{Z}_j(\omega ),$$

$$x_{k,j}={\mu }_{k,j}+\frac{2\pi }{N}\sum _{n=1}^{N-1}(c_{k,j,n}\cos ({\omega }_{n}k)-d_{k,j,n}\sin ({\omega }_{n}k)),$$

$${\psi }_{k,j}({\omega }_n)d{\omega }_n=|A_{k,j}({\omega }_n){|}^2\mathbb{E}|d{Z}_j({\omega }_n){|}^2.$$

$${\psi }_{k,j}({\omega }_n)=\frac{\pi }{N}\mathbb{E}[(c_{k,j,n}+i{d}_{k,j,n})(c_{k,j,n}-i{d}_{k,j,n})].$$

$${\mathbf{\Psi }}_k({\omega }_n)=\frac{\pi }{N}\mathbb{E}[({\mathbf{c}}_{k,n}+i{\mathbf{d}}_{k,n})({\mathbf{c}}_{k,n}-i{\mathbf{d}}_{k,n}{)}^{\top }],$$

$$x_{k,j}={\mu }_{m,j}+\frac{2\pi }{N}\sum _{n=1}^{N-1}(p_{m,j,n}\cos ({\omega }_{n}k)-q_{m,j,n}\sin ({\omega }_{n}k)),$$

$$\displaystyle\mathbf{A}_{m}:=\frac{2\pi}{N}\leavevmode\resizebox{368.57964pt}{}{$\begin{bmatrix}1&\cos\left(\frac{(m-1)W+1}{N}\pi\right)&-\sin\left(\frac{(m-1)W+1}{N}\pi\right)&\dots&\cos\left(\frac{(N-1)((m-1)W+1)}{N}\pi\right)&-\sin\left(\frac{(N-1)((m-1)W+1)}{N}\pi\right)\\ 1&\cos\left(\frac{(m-1)W+2}{N}\pi\right)&-\sin\left(\frac{(m-1)W+2}{N}\pi\right)&\dots&\cos\left(\frac{(N-1)((m-1)W+2)}{N}\pi\right)&-\sin\left(\frac{(N-1)((m-1)W+2)}{N}\pi\right)\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 1&\cos\left(\frac{mW}{N}\pi\right)&-\sin\left(\frac{mW}{N}\pi\right)&\dots&\cos\left(\frac{(N-1)mW}{N}\pi\right)&-\sin\left(\frac{(N-1)mW}{N}\pi\right)\\ \end{bmatrix}$}.$$

$${\mathbf{\Psi }}_m({\omega }_n)=\frac{\pi }{N}\mathbb{E}[({\mathbf{p}}_{m,n}+i{\mathbf{q}}_{m,n})({\mathbf{p}}_{m,n}-i{\mathbf{q}}_{m,n}{)}^{\top }],$$

$$\begin{array}{rl}{S}^{\mathsf{mt}}(\omega )& =\frac{1}{P}\sum _{p=1}^P|y^{(p)}(\omega ){|}^2\text{with},\\ y^{(p)}(\omega )& =\sum _{k=1}^K{\nu }_k^{(p)}{x}_k{e}^{-i\omega k},\end{array}$$

$$\begin{array}{rl}{S}_{r,t}^{\mathsf{mt}}(\omega )& =\frac{1}{P}\sum _{p=1}^P{y}_r^{(p)}(\omega )(y_t^{(p)}(\omega ){)}^{\ast }\text{where,}\\ y_r^{(p)}(\omega )& =\sum _{k=1}^K{\nu }_k^{(p)}{x}_{k,r}{e}^{-i\omega k}.\end{array}$$

$$\log f(\mathcal{D}|\mathbf{X})$$

$$\displaystyle=\sum_{j=1}^{J}\sum_{k=1}^{K}L\Big{\{}\overline{n}_{k,j}x_{k,j}\!-\!\log(1\!+\!\exp(x_{k,j})\Big{\}}$$

$$\begin{array}{r}{\overline{n}}_{k,j}=\frac{1}{L}\sum _{l=1}^L{n}_{k,j}^{(l)}\to \mathrm{logistic}(x_{k,j})a.s.\end{array}$$

$$\begin{array}{r}{\frac{\partial \log f(\{n_{k,j}^{(l)}{\}}_{l=1}^L|x_{k,j})}{\partial x_{k,j}}|}_{x_{k,j}=\mathrm{logit}({\overline{n}}_{k,j})}=0,\end{array}$$

$$\begin{array}{r}({\overline{n}}_{k,j}{)}^{(p)}=\mathrm{logistic}({\nu }_{k\mathrm{mod}W+1}^{(p)}\mathrm{logit}({\overline{n}}_{k,j})),\end{array}$$

$${\mathbf{w}}_m=\mathbf{\Phi }{\mathbf{w}}_{m-1}+{\mathbf{\eta }}_m$$

$$\begin{split}&\log\,f(\mathcal{D},\mathbf{V}|\bm{\theta})=\log\,f(\mathcal{D}|\mathbf{V},\bm{\theta})\>+\>\log\,f_{\mathbf{V}|\bm{\theta}}(\mathbf{V}|\bm{\theta})\\ &=\sum_{j=1}^{J}\sum_{m=1}^{M}\sum_{w=1}^{W}\sum_{l=1}^{L}\Big{\{}n_{(m-1)W+w,j}^{(l)}(\mathbf{A}_{m}\mathbf{V}_{m})_{w,j}\\ &-\log\,(1+\exp(\mathbf{A}_{m}\mathbf{V}_{m})_{w,j}))\Big{\}}-\frac{1}{2}\sum_{m=1}^{M}\,\Big{\{}\log|\mathbf{Q}_{m}|\\ &+(\mathbf{w}_{m}-\mathbf{\Phi}\mathbf{w}_{m-1})^{\top}\mathbf{Q}_{m}^{-1}(\mathbf{w}_{m}-\mathbf{\Phi}\mathbf{w}_{m-1})\Big{\}}+C_{1},\end{split}$$

$$\begin{split}f(\mathbf{Q}_{m})\propto\leavevmode\resizebox{177.78578pt}{}{$\displaystyle\exp\!\left(\!\!-\rho\sum_{j=1}^{2J}\sum_{n=1}^{N-2}\Big{(}\log(Q_{m,J(2n-1)+j})\!-\!\log(Q_{m,J(2n+1)+j})\Big{)}^{2}\!\right)$}.\end{split}$$

$$\log$$

$$\displaystyle\leavevmode\resizebox{203.80193pt}{}{$\displaystyle-\rho\sum_{m=1}^{M}\sum_{j=1}^{2J}\sum_{n=1}^{N-2}\Big{(}\log(Q_{m,J(2n-1)+j})-\log(Q_{m,J(2n+1)+j})\Big{)}^{2}+C_{3}$},$$

$${\mathcal{Q}}^{(r)}:=\mathbb{E}[\log f(\mathcal{D},\mathbf{V},\mathbf{\theta })|\mathcal{D},{\hat{\mathbf{\theta }}}^{(r)}]$$

$${\mathbf{w}}_{m|M}:=\mathbb{E}[{\mathbf{w}}_m|\mathcal{D},{\hat{\mathbf{\theta }}}^{(r)}]$$

$${\mathbf{\Sigma }}_{m|M}:=\mathbb{E}[({\mathbf{w}}_m-{\mathbf{w}}_{m|M})({\mathbf{w}}_m-{\mathbf{w}}_{m|M}{)}^{\top }|\mathcal{D},{\hat{\mathbf{\theta }}}^{(r)}]$$

$${\mathbf{\Sigma }}_{m,m-1|M}:=\mathbb{E}[({\mathbf{w}}_m-{\mathbf{w}}_{m|M})({\mathbf{w}}_{m-1}-{\mathbf{w}}_{m-1|M}{)}^{\top }|\mathcal{D},{\hat{\mathbf{\theta }}}^{(r)}]$$

$${\mathbf{w}}_{m|m}^{(r)}$$

$$\displaystyle=\underset{\mathbf{w}_{m}}{\textrm{argmax}}\bigg{(}\sum_{j=1}^{J}\sum_{s=1}^{m}\sum_{w=1}^{W}L\,\Big{\{}\overline{n}_{(s-1)W+w,j}(\mathbf{A}_{s}\mathbf{V}_{s})_{w,j}$$

$$\displaystyle-\log\,(1+\exp(\mathbf{A}_{s}\mathbf{V}_{s})_{w,j}))\big{\}}-\frac{1}{2}\sum_{s=1}^{m}\,\Big{\{}\log|\mathbf{Q}_{s}^{(r)}|$$

$$\displaystyle+(\mathbf{w}_{s}\!-\!\mathbf{\Phi}\mathbf{w}_{s-1})^{\top}(\mathbf{Q}_{s}^{(r)})^{-1}(\mathbf{w}_{s}\!-\!\mathbf{\Phi}\mathbf{w}_{s-1})\Big{\}}\bigg{)},$$

$$\begin{array}{rl}& {\mathbf{Q}}_m^{(r+1)}=\underset{{\mathbf{Q}}_m}{\mathrm{argmax}}{\mathcal{Q}}^{(r)}.\end{array}$$

$$\frac{\partial {\mathcal{Q}}^{(r)}}{\partial Q_{m,i}}=-\frac{1}{2}\left\{\frac{1}{Q_{m,i}}-\frac{P_{m,i}}{Q_{m,i}^2}\right\}-\frac{2\rho }{Q_{m,i}}\log \left(\frac{Q_{m,i}^2}{Q_{m,i-2J}{Q}_{m,i+2J}}\right),$$

$${\mathbf{P}}_m:=$$

$$=$$

Full text

I Introduction

Neural oscillations are known to play a significant role in mediating the cognitive and motor functions of the brain [2, 3, 4]. The advent of high-density electrophysiology recordings [5, 6, 7] from multiple locations in the brain has opened a unique window of opportunity to probe these oscillations at the neuronal scale. In order to exploit such experimental data for inferring the mechanisms of brain function, spectral analysis techniques tailored for such neuronal spiking data are required [8].

Existing techniques for spectral analysis of neuronal data use point process theory [8, 9, 10] to capture the spiking statistics. Due to the time-domain smoothing procedures used by existing techniques [11, 12, 13] for recovering the latent processes that drive spiking activity, the resulting power spectral density (PSD) estimates undergo distortion in the spectral domain. Alternative methods to directly estimate the PSD from spiking data have recently been proposed in [14, 15].

These existing methods consider univariate spiking observations and assume the latent process to be second-order stationary during the observation period. However, it is known that the brain oscillations that drive neuronal spiking are non-stationary and may exhibit rapid changes corresponding to the brain state or behavioral dynamics [12, 16].

Non-stationary time series analysis has been well studied for multivariate continuous signals and various methods have been proposed to quantify the energy-frequency-time distributions [17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. One notable example is the evolutionary power spectral characterization [22, 23], which defines a non-stationary spectral density matrix in order to quantify the local spectral energy distributions at each instant of time for a multivariate process. This characterization has a physical interpretation similar to that of stationary Fourier spectral analysis, and reduces to the classical power spectra when the processes are stationary [22]. A unified approach that considers multivariate spiking observations driven by non-stationary latent processes is lacking, but highly desired due to the emerging demands of modern neuronal data analysis.

In this work, we close this gap by developing a framework to estimate the evolutionary spectral density (ESD) matrix of a multivariate non-stationary latent process, given spiking observations. We model the spiking observations as multiple realizations of point processes with logistic links to the latent continuous processes. We then pose the problem of spectral estimation within a multitapering framework. Multitapering is a widely-used PSD estimation technique with desirable bias-variance trade-off performance [27], due to the usage of the discrete prolate spheroidal sequences (dpss) [28] as data tapers, which are known for their minimal spectral leakage [29, 30, 31].

We employ a state-space model to characterize the dynamics of the evolutionary spectra, with the underlying states pertaining to the eigen-spectra of the multivariate latent processes [25]. We derive an Expectation-Maximization (EM) algorithm for efficiently computing the maximum a posteriori (MAP) estimate of the latent variables and smoothed states given the spiking observations, which we then use to construct the ESD matrix. We establish theoretical bounds on the bias-variance performance of our proposed estimator, and recover the favorable asymptotic properties of the classical multitaper framework.

We compare the performance of our proposed method to existing techniques through simulation studies and application to experimentally-recorded neuronal data. We present two simulated case studies using non-stationary multivariate autoregressive processes, whose dynamics are inspired by neural oscillations. These studies reveal that the proposed method outperforms two of the widely used methods for deriving spectral representations from spiking data. Finally, we apply our proposed estimator to multi-unit spike and local field potential (LFP) recordings from a human subject undergoing general anesthesia [12, 14]. Our proposed method corroborates existing hypotheses on the relation between the LFP signals and spiking dynamics, by providing a high resolution characterization of the underlying spectrotemporal couplings.

The rest of the paper is organized as follows: We present our problem formulation in Section II, followed by the proposed estimation framework in Section III. We provide our theoretical results on the bias-variance performance of the proposed estimator in Section IV. Our simulation studies are presented in Section V, followed by application to experimentally-recorded data in Section VI. Finally, we close the paper by our concluding remarks in Section VII.

II Problem Formulation

Let $N(t)$ and $H(t)$ denote the point process representing the number of spikes and spiking history of a neuron in $[0,t)$, respectively, where $t\in[0,T]$ and $T$ denotes the observation duration. The Conditional Intensity Function (CIF) [32] 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$ is an integer (with no loss of generality). We refer to $\lambda_{k}$ as CIF hereafter for brevity.

In a similar fashion, we consider spiking observations from an ensemble of $J$ neurons, with CIFs $\{\lambda_{k,j}\}_{k=1}^{K}$, for $j=1,2,\cdots,J$. Suppose that for each neuron, $L$ independent realizations of the spiking activity are observed. The collection of the binary spiking observations are represented as $\{n_{k,j}^{(l)}\}_{l,k,j=1}^{L,K,J}$. We model the $j^{\sf th}$ CIF by a logistic link to a latent random process, $\mathbf{x}_{j}=[x_{1,j},x_{2,j},\cdots,x_{K,j}]^{\top}$, which needs not be stationary in general. Thus, the matrix $\mathbf{X}=[\mathbf{x}_{1},\;\mathbf{x}_{2},\cdots,\mathbf{x}_{J}]_{K\times J}$ represents a $J$-variate random process.

Accordingly, for $1\leq j\leq J,\;1\leq k\leq K$ and $1\leq l\leq L$, we have

[TABLE]

where $\lambda_{k,j}=\operatorname{logistic}(x_{k,j})=1/(1+\exp{(-x_{k,j})})$. Our goal is to estimate the time variant power spectral density of each process $\mathbf{x}_{j}$, for $1\leq j\leq J$, and the time variant cross spectra between each pair of processes. Following the formation of Priestley’s Evolutionary Spectra [22], each random process $x_{k,j}$, with mean $\mu_{k,j}$, will have a representation of the form,

[TABLE]

where $A_{k,j}(\omega)$ is the time-varying amplitude function and $dZ_{j}(\omega)$ is an orthogonal increment process. To define a discrete-parameter harmonic process, we approximate $Z_{j}(\omega)$ by a jump process over $N$ frequency bins [14], and thereby replace $dZ_{j}(\omega_{n})$ with $\frac{\pi}{N}(a_{j,n}+ib_{j,n})$, at $\omega_{n}=n\pi/N,1\leq n\leq N-1$. Given the random process is real, using the symmetry $Z_{j}(\omega)=Z_{j}(-\omega)$, we express the discretization of Eq. (3) as

[TABLE]

where $c_{k,j,n}=A_{k,j}(\omega_{n})\>a_{j,n}$ and $d_{k,j,n}=A_{k,j}(\omega_{n})\>b_{j,n}$ are real-valued random variables. Further, the evolutionary spectrum [22] at time $k$ is defined as

[TABLE]

Hence, according to our model, the ESD function can be approximated by,

[TABLE]

The Spectral Density Matrix of a multivariate stationary random process is defined as, $\bm{\Psi}(\omega)=\frac{1}{2\pi}\sum_{\ell=-\infty}^{\infty}e^{-i\ell\omega}\,\bm{\Gamma}(\ell)$, where $\bm{\Gamma}(\cdot)$ is the covariance matrix of the process [33]. Further, considering a vector valued orthogonal increment process $\mathbf{Z}(\omega)=[Z_{1}(\omega),Z_{2}(\omega),\cdots,Z_{J}(\omega)]^{\top}$, the spectral density matrix is characterized by $\bm{\Psi}(\omega)d\omega=\mathbb{E}[d\mathbf{Z}(\omega)d\mathbf{Z}(\omega)^{H}]$. We extend this to the evolutionary spectra, and formulate the ESD matrix at time $k$ and frequency $\omega_{n}$ for our model as,

[TABLE]

where $\mathbf{c}_{k,n}=[c_{k,1,n},\,c_{k,2,n},\cdots,c_{k,J,n}]^{\top}$ and $\mathbf{d}_{k,n}=[d_{k,1,n},\,d_{k,2,n},\cdots,d_{k,J,n}]^{\top}$.

Further, we assume the processes to be quasi-stationary [25], the $J$-variate random process to be jointly stationary in windows of length $W$. The total data duration, $K$ is divided into $M$ non overlapping segments of length $W$, with $K=MW$. Thus, the vector process $[x_{k,1},\,x_{k,2},\cdots,x_{k,J}]$ is assumed to be jointly stationary for $(m-1)W+1\leq k\leq mW$, $1\leq m\leq M$. Simplifying the previous model under this quasi-stationarity assumption, for $(m-1)W+1\leq k\leq mW$, $1\leq m\leq M$ we get,

[TABLE]

where $p_{m,j,n}$ and $q_{m,j,n}$ are real-valued random variables.

Defining $\mathbf{X}_{m,j}=[x_{(m-1)W+1,j},x_{(m-1)W+2,j},\cdots,$ $x_{mW,j}]^{\top}$, $\mathbf{v}_{m,j}=[\frac{N}{2\pi}\mu_{m,j},\;p_{m,j,1},\;q_{m,j,1},\cdots,p_{m,j,N-1},$ $q_{m,j,N-1}]^{\top}$, $\bar{\mathbf{X}}_{m}=[\mathbf{X}_{m,1},\mathbf{X}_{m,2},\cdots,\mathbf{X}_{m,J}]_{W\times J}$, $\mathbf{V}_{m}=[\mathbf{v}_{m,1},\mathbf{v}_{m,2},\cdots,\mathbf{v}_{m,J}]_{(2N-1)\times J}$ and $\mathbf{A}_{m}$ as in Eq. (4), we get $\mathbf{X}_{m,j}=\mathbf{A}_{m}\mathbf{v}_{m,j}$ and $\bar{\mathbf{X}}_{m}=\mathbf{A}_{m}\mathbf{V}_{m}$. Further, we define $\mathbf{p}_{m,n}=[p_{m,1,n},p_{m,2,n},\cdots,p_{m,J,n}]^{\top}$, $\mathbf{q}_{m,n}=[q_{m,1,n},q_{m,2,n},\cdots,q_{m,J,n}]^{\top}$, $\mathbf{w}_{m,0}=[\mu_{m,1},\mu_{m,2},\cdots,\mu_{m,J}]^{\top}$, $\mathbf{w}_{m,n}=[\mathbf{p}_{m,n}^{\top}\,,\,\mathbf{q}_{m,n}^{\top}]^{\top}_{2J\times 1}$ for $1\leq n\leq N-1$ and $\mathbf{w}_{m}=[\mathbf{w}_{m,0}^{\top},\mathbf{w}_{m,1}^{\top},\cdots,\mathbf{w}_{m,N-1}^{\top}]^{\top}_{J(2N-1)\times 1}$. Note that $\mathbf{w}_{m}$ is the vectorization of the matrix $\mathbf{V}_{m}$ and both are equivalent representations, for the discrete parameter harmonic process that drives the spiking observations in the time window $m$.

Accordingly, the ESD matrix,

[TABLE]

can be determined if $\mathbb{E}[\mathbf{w}_{m,n}\mathbf{w}_{m,n}^{\top}]$ is estimated for $1\leq n\leq N-1$. Thus, the task of determining the evolutionary power spectra of the $J$-variate random process can be reduced to computing $\mathbb{E}[\mathbf{w}_{m}\mathbf{w}_{m}^{\top}]$, for $m=1,2,\cdots,M$, given the ensemble of spiking data $\mathcal{D}=\{n_{k,j}^{(l)}\}_{l,k,j=1}^{L,K,J}$.

III Multitaper Estimate of the Spectral Density Matrix

III-A The Multitaper Framework

Multitapering is a technique widely used in power spectral estimation of stationary random processes, to overcome bias and variance limitations of conventional Fourier analysis [27]. The multitaper spectral estimate of a stationary process, $x_{1},x_{2},\cdots,x_{K}$ is defined as,

[TABLE]

where $\nu_{k}^{(p)}$ is the $k^{\sf th}$ time sample of the $p^{\sf th}$ discrete prolate spheroidal sequence (dpss) [28], for $1\leq p\leq P$. The dpss tapers are a set of orthogonal tapers that maximally concentrate their spectral power within a design bandwidth of $[-\frac{\xi f_{s}}{K},\frac{\xi f_{s}}{K}]$, for some positive constant $\xi$. While originally designed for univariate time series, multitapering has been likewise extended to multivariate time series [34]. Given a second order jointly stationary $J$-dimensional random process $\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{K}$, where, $\mathbf{x}_{k}=[x_{k,1},x_{k,2},\cdots,x_{k,J}]$, the multitaper cross-spectral estimate between the $r^{\sf th}$ process and the $t^{\sf th}$ process $(r,t\in\{1,2,\cdots,J\})$ has been defined as,

[TABLE]

The spectral estimate $y_{r}^{(p)}(\omega)(y_{t}^{(p)}(\omega))^{*}$ corresponding to the $p^{\sf th}$ taper is referred to as the $p^{\sf th}$ eigen-spectral estimate. Given that our observed data are binary random variables parameterized by a function of $\mathbf{X}$ (and not the actual process $\mathbf{X}$), determining the multitaper spectral estimate is not straightforward. However, note that the data log-likelihood,

[TABLE]

depends on the observations only through the ensemble mean, $\overline{n}_{k,j}:=\frac{1}{L}\sum_{l=1}^{L}n_{k,j}^{(l)}$, resulting $\overline{n}_{k,j}$ to be a sufficient statistic. Thus, if one could characterize the effect of tapering the time series on the ensemble mean, it would be possible to form the multitaper spectral estimate of the unobserved time series.

Given that $n_{k,j}^{(l)}\sim\operatorname{Bernoulli}(\lambda_{k,j})$, with $\lambda_{k,j}=\operatorname{logistic}(x_{k,j}):=1/(1+\exp(-x_{k,j}))$, we have $x_{k,j}=\operatorname{logit}(\lambda_{k,j}):=\log(\lambda_{k,j}/(1-\lambda_{k,j}))$ and,

[TABLE]

by the strong law of large numbers, for $\overline{n}_{k,j}\neq 0,1$. Further, consider the log-likelihood \log f(\{n_{k,j}^{(l)}\}_{l=1}^{L}|x_{k,j})=L\Big{\{}\overline{n}_{k,j}x_{k,j}-\log\,(1+\exp(x_{k,j})\Big{\}}, in Eq. (8). We notice that

[TABLE]

which implies that $\operatorname{logit}(\overline{n}_{k,j})$ is the maximum likelihood estimator of $x_{k,j}$. Thus, we take $\operatorname{logit}(\overline{n}_{k,j})$ as an estimator of $x_{k,j}$. By extending this argument to the tapered processes $\{\nu_{k\!\!\!\mod\!W+1}^{(p)}x_{k,j}\}_{k,j,p=1}^{K,J,P}$, we can similarly take

[TABLE]

to be the estimator of the ensemble mean that would have been generated if the random process had been filtered by the $p^{\sf th}$ data taper, for $1\leq k\leq K$. Note that when $\overline{n}_{k,j}=0$ or $\overline{n}_{k,j}=1$, the function $\operatorname{logit}(\overline{n}_{k,j})$ is not defined. In such cases, we directly estimate $(\overline{n}_{k,j})^{(p)}$ by $\overline{n}_{k,j}$.

We thus need to compute the evolutionary spectra corresponding to each of the $P$ tapers, and finally derive the multitaper estimate by averaging the $P$ eigen-spectral estimates. In the next subsection, we consider estimating the power spectral density matrix of the untapered process first, and then extend it to the case of $P$ tapers by replacing the ensemble average of spiking data $\{\overline{n}_{k,j}\}$ with the tapered ensemble mean $\{\overline{n}_{k,j}^{(p)}\}$, for $p=1,2,\cdots,P$.

III-B MAP Estimation of the Parameters via the EM Algorithm

In order to capture the evolution of the spectra, we impose a stochastic continuity constraint on the random variables $\mathbf{w}_{m}$, following the discrete state-space model,

[TABLE]

where the state transition matrix $\mathbf{\Phi}$ is a constant matrix and $\bm{\eta}_{m}\sim\mathcal{N}(\mathbf{0},\mathbf{Q}_{m})$. We consider the special case where $\mathbf{\Phi}$ takes the form $\alpha\mathbf{I}$, for simplicity. Nevertheless, $\mathbf{\Phi}$ can also be estimated from data within the same Expectation-Maximization framework that follows next [35].

The parameters to be estimated are $\bm{\theta}:=\{\mathbf{Q}_{m},\,1\leq m\leq M\}$, and the observations are binary spiking data $\mathcal{D}=\{n_{k,j}^{(l)}\}_{k,j,l\,=1}^{K,J,L}$. Considering the variable $\mathbf{V}=\{\mathbf{V}_{m},\,1\leq m\leq M\}$ to be the hidden data, we aim at recovering $\bm{\theta}$ via a MAP estimation problem. First, the log-likelihood of the complete data has the form

[TABLE]

where $\mathbf{w}_{0}=\mathbf{0}$ and $C_{1}$ represents terms that are not functions of $\mathcal{D},\mathbf{V}$ or $\bm{\theta}$. Next, considering that neuronal spiking activity is typically sparse in time and the number of observed realizations ($L$) is limited, an appropriate prior distribution on the parameters helps in reducing the estimation variance. We consider a diagonal covariance matrix $\mathbf{Q}_{m}$, whose $i^{\sf th}$ diagonal entry is denoted by $Q_{m,i}$. Further, we assume $\mathbf{Q}_{m}$ to be independent and identically distributed for $1\leq m\leq M$, with a distribution of the form,

[TABLE]

This prior distribution encourages continuity in log scale of the spectral estimates corresponding to adjacent frequency bins of each CIF, and can be controlled by appropriately selecting the hyper-parameter $\rho$. Accordingly, considering $(\mathcal{D},\mathbf{V})$ to be the set of complete data, the joint distribution of the MAP estimation problem is formulated as,

[TABLE]

where $C_{2}$ and $C_{3}$ represent terms that do not depend on $\bm{\theta}$. We next construct the EM algorithm for solving the MAP estimation problem:

III-B1 E Step

Suppose that the current estimate of $\bm{\theta}$ at the $r^{\sf th}$ iteration is given by $\widehat{\bm{\theta}}^{(r)}$. Then, the Q-function

[TABLE]

can be evaluated if the conditional expectations

[TABLE]

are known. To compute the conditional expectations $\mathbf{w}_{m|M}$, $\mathbf{\Sigma}_{m|M}$ and $\mathbf{\Sigma}_{m,m-1|M}$ we utilize the Fixed Interval Smoothing [36] and the Covariance Smoothing [37] algorithms. However, considering that the forward model is not Gaussian in our formulation, we cannot directly use these algorithms to estimate $\mathbf{w}_{m|m}$ and $\mathbf{\Sigma}_{m|m}$.

Hence, we employ an alternative method to estimate these conditional moments, utilizing the distribution $f({\{\mathbf{V}_{s}\}_{s=1}^{m}|\mathcal{D}_{1}^{m},\hat{\bm{\theta}}^{(r)}})$, where $\mathcal{D}_{1}^{m}=\{n_{k,j}^{(l)}\}_{k,j,l\,=1}^{mW,J,L}$. From Bayes’s theorem we see that this is proportional to the product of the two distributions, $f(\mathcal{D}_{1}^{m}|\{\mathbf{V}_{s}\}_{s=1}^{m},\hat{\bm{\theta}}^{(r)})$ and $f(\{\mathbf{V}_{s}\}_{s=1}^{m}|\hat{\bm{\theta}}^{(r)})$, which are Binomial and Gaussian distributed, respectively. Accordingly, we see that the distribution $\{\mathbf{V}_{s}\}_{s=1}^{m}|\mathcal{D}_{1}^{m},\hat{\bm{\theta}}^{(r)}$ is unimodal, and hence we approximate it by a multivariate Gaussian, and derive the mean of the distribution, $\mathbf{w}_{m|m}^{(r)}$ by the mode of $\log\,f({\{\mathbf{V}_{s}\}_{s=1}^{m}|\mathcal{D}_{1}^{m},\hat{\bm{\theta}}^{(r)}})$:

[TABLE]

and its covariance, $\mathbf{\Sigma}_{m|m}^{(r)}$ by the negative of the inverse Hessian of $\log f({\{\mathbf{V}_{s}\}_{s=1}^{m}|\mathcal{D}_{1}^{m},\hat{\bm{\theta}}^{(r)}})$. Observing that the objective function is a combination of convex functions and is differentiable, we perform the above optimization for $\mathbf{w}_{m|m}^{(r)}$ using the Newton-Raphson method. Further, we concurrently estimate $\mathbf{\Sigma}_{m|m}^{(r)}$ using the Hessian matrix of the objective.

III-B2 M Step

Noticing that the Q-function $\mathcal{Q}^{(r)}$ is separable in terms of $\mathbf{Q}_{m}$’s, we can update $\mathbf{Q}_{m}$ independently, for $1\leq m\leq M$:

[TABLE]

The partial derivative of the log-likelihood function with respect to each diagonal element of $\mathbf{Q}_{m}$ takes the form,

[TABLE]

where $P_{m,i}$ is the $i^{\sf th}$ diagonal element of the matrix

[TABLE]

Accordingly, we employ the multivariate Newton-Raphson algorithm to perform this maximization and derive the updates for $\mathbf{Q}_{m}$, $1\leq m\leq M$. Following convergence, we use the final estimates of $\mathbf{w}_{m|M}$ and $\mathbf{\Sigma}_{m|M}$ derived through the above EM procedure to compute

[TABLE]

for $1\leq m\leq M$. Then, for $1\leq n\leq N-1$, the ESD matrix according to Eq. (5) is estimated as

[TABLE]

where $\mathbf{R}_{m}^{n}:=\mathbb{E}[\mathbf{w}_{m,n}\mathbf{w}_{m,n}^{\top}|\mathcal{D},\bm{\theta}]$ is a submatrix of $\mathbf{R}_{m}$ defined as:

[TABLE]

An implementation of this estimation procedure is outlined in Algorithm 1. As mentioned in Section III-A, the same EM procedure can be carried out for $\{\overline{n}_{k,j}^{(p)}\}$, for $p=1,2,\cdots,P$ in order to estimate the multivariate eigen-spectra. Finally, the multitaper spectral density matrix is formed by averaging the eigen-spectral estimates as outlined in Algorithm 2. We refer to our proposed algorithm as the Point Process Multitaper Evolutionary Spectral Density (PPMT-ESD) estimator.

IV Theoretical Analysis

In this section we derive bounds on the bias and variance of the proposed PPMT-ESD estimator. We first briefly review the corresponding bounds for the direct multitaper estimator (Eq. (6)) of a time-series. As proven in [38], the bias and variance of the multitaper estimate of a stationary process $x_{1},x_{2},\cdots,x_{K}$ with a uniformly continuous PSD $S(\omega)$ are bounded as follows:

[TABLE]

as $K\rightarrow\infty$. Here $c_{p}$ is the eigenvalue corresponding to the taper $\nu^{(p)}$ and $\beta(\omega)=0$ if $\frac{\omega}{2\pi}\neq 0,1/2\,\operatorname{mod}1$ and is equal to $1$ otherwise. It is evident that the multitaper estimator $\widehat{S}^{\sf mt}(\omega)$ is asymptotically unbiased.

We state our main theorem for a univariate second-order stationary process $x_{1},x_{2},\cdots,x_{K}$, corresponding to the special case of $J=1$, for the clarity of exposition. We later on provide extensions to the multivariate and quasi-stationary cases. In order to apply the same treatment to the case of multitaper estimate of the evolutionary spectra from spiking observations, we need to make two extra technical assumptions.

Assumption (1). From Eq. (8), the data likelihood for the univariate case can be expressed as

[TABLE]

where $\overline{n}_{k}:=\frac{1}{L}\sum_{l=1}^{L}n_{k}^{(l)}$. Given the nonlinear functional form of the integrand, we consider the saddle point approximation [39] of the integral, and take $\operatorname{logit}(\overline{n}_{k})$ as an estimator of $x_{k}$. Under this approximation, the multitaper spectral estimate is given by

[TABLE]

[TABLE]

Assumption (2). We assume that $|x_{k}|\;\leq B$ for all $k=1,2,\cdots,K$, for some fixed upper bound $B$. In defining the bias and variance, we condition the expectations on the event $A:=\{\overline{n}_{k}\;|\;\overline{n}_{k}\neq 0,\,\overline{n}_{k}\neq 1,\,1\leq k\leq K\}$, which is highly probable due to the absolute boundedness of $x_{k}$, for large $L$ (See Appendix A for details). We denote the conditional bias and variance by $\operatorname{bias}_{A}(\cdot)$ and $\operatorname{Var}_{A}(\cdot)$, respectively. Note that for the multivariate case, we naturally extend this assumption to $|x_{k,j}|\leq B$ for all $k,j$, and define the set $A$ as $A:=\{\overline{n}_{k,j}\;|\;\overline{n}_{k,j}\neq 0,\,\overline{n}_{k,j}\neq 1,\,1\leq k\leq K,1\leq j\leq J\}$.

Theorem 1** (Univariate Case).**

Under the Assumptions (1) and (2) and for sufficiently large $L$, the conditional bias and variance of $\widehat{S}^{\sf mt}(\omega)$ in Eq. (19) is bounded with respect to those of the direct multitaper estimate $S^{\sf mt}(\omega)$ given in Eq. (18) as:

[TABLE]

[TABLE]

where $g_{1}$ and $g_{2}$ are bounded constants depending on $B$, $K$ and $L$.

Proof of Theorem 1.

The proof of Theorem 1 is given in Appendix A. ∎

Remark. In words, Theorem 1 states that the cost of the indirect access to the process $\{x_{k}\}_{k=1}^{K}$ through spiking observations with $L$ trials appears as excess terms in both the bias and variance, which would go to zero as $\frac{L}{K^{2}\log^{2}L}\rightarrow\infty$. Hence, for large enough number of realizations $L$, one can expect a performance close to the direct multitaper estimate of $\{x_{k}\}_{k=1}^{K}$. The result of Theorem 1 can be extended to the case of unconditional expectations, if the inverse mapping $\max(B,\,\min(-B,\,\operatorname{logit}(\overline{n}_{k})))$ is adopted instead of $\operatorname{logit}(\overline{n}_{k})$. Then, the bounds follow directly even without conditioning on $A$. While Assumption (2) on the boundedness of the time-series is natural in practical scenarios, as long as $B\leq\epsilon\log L$, for some $\epsilon<1/2$, the excess bias and variance terms will be bounded by $\mathcal{O}(\log^{2}L/L^{1/2-\epsilon})$, which will converge to zero asymptotically. Thus, the bias and variance of the estimator $\widehat{S}^{\sf mt}(\omega)$ will converge to those of $S^{\sf mt}(\omega)$ even under the milder condition, $|x_{k}|\;\leq\,\epsilon\log L$ for $1\leq k\leq K$, and some $\epsilon<1/2$.

The following corollary extends Theorem 1 to the multivariate case:

Corollary 1** (Stationary Multivariate Case).**

Consider a second order jointly stationary J-variate process $\{\mathbf{x}_{k}\}_{k=1}^{K}$, where, $\mathbf{x}_{k}=[x_{k,1},x_{k,2},\cdots,x_{k,J}]^{\top}$. Suppose that the observations are binary spiking data, $\{n_{k,j}^{(l)}\}_{k,j,l=1}^{K,J,L}$ with $n_{k,j}^{(l)}\sim\operatorname{Bernoulli}({\lambda}_{k,j})$ and $\overline{n}_{k,j}:=\frac{1}{L}\sum_{l=1}^{L}n_{k,j}^{(l)}$, where ${\lambda}_{k,j}=\operatorname{logistic}\,(x_{k,j})$. Then, under Assumptions (1) and (2), the bias and variance of the multitaper cross-spectral estimate between the $r^{\sf th}$ and $t^{\sf th}$ processes, given by

[TABLE]

are bounded as follows:

[TABLE]

[TABLE]

where $g^{\prime}_{1}$ and $g^{\prime}_{2}$ are bounded constants and depend only on $B$, $K$, and $L$.

Proof.

The proof is given in Appendix B. ∎

It is not difficult to verify that $\operatorname{bias}(S^{\sf mt}_{r,t}(\omega))$ and $\operatorname{Var}(S^{\sf mt}_{r,t}(\omega))$ can be upper bounded in a similar fashion to Eq. (18), with the true cross-spectra $S_{r,t}(\omega)$ replacing $S(\omega)$. Before extending the result of Corollary 1 to the quasi-stationary case, we need an additional assumption:

Assumption (3). Given that Corollary 1 holds for large $L$, in this regime we relax the prior distribution on $\mathbf{Q}_{m}$ to be flat, i.e., $f(\mathbf{Q}_{m})\propto 1$. Recall that the rationale for using a prior on $\mathbf{Q}_{m}$ in Section III-B was to reduce the variance of the estimates in the low spiking regime, i.e., small $L$.

Finally, combining Corollary 1 and the treatment of [25], we have the following corollary on the bias and variance of the PPMT-ESD estimator:

Corollary 2** (Quasi-stationary Multivariate Case).**

Suppose that the J-variate process in Corollary 1 is quasi-stationary (jointly stationary within consecutive windows of length $W$). Let $\Psi_{m,r,t}(\omega_{n})$ be the cross-spectra between the $r^{\sf th}$ and $t^{\sf th}$ processes over window $m$, $1\leq m\leq K/W$, and $\widehat{\Psi}^{\sf mt}_{m,r,t}(\omega_{n})$ be the corresponding multitaper estimate obtained from spiking observations. Then, under Assumptions (1)–(3), the bias and variance of the proposed PPMT-ESD estimator at window $m$ can be bounded as,

[TABLE]

[TABLE]

where $g^{\prime\prime}_{1}(\omega)$, $g^{\prime\prime}_{2}(\omega)$ are bounded functions of $B$, $L$, $W$ and $\kappa_{m}(\cdot)$ is a function of $\omega$, given in Appendix B.

Proof.

The proof is given in Appendix B. ∎

V Simulation Studies

First, we applying the PPMT-ESD estimator to simulated data and compare its performance to some of the existing methods. We consider two cases based on multivariate autoregressive processes:

Case 1: Estimating the spectral density matrix of a latent trivariate random process, given binary spiking observations.

Case 2: Estimating the spectral density matrix of two processes, where one is directly observable and the other is latent and observed via binary spiking.

As for comparison, we consider two existing methods for extracting spectral representations of spiking data:

State-Space ESD Estimator

This estimator closely follows the approach in [11]. Recall that the observations are $n_{k}^{(l)}\sim\operatorname{Bernoulli}({\lambda}_{k})$, where ${\lambda}_{k}=\operatorname{logistic}\,(x_{k})$. Using the same notation as in [11], we model the latent process $x_{k}$ as a first-order autoregressive model, $x_{k}=x_{k-1}+\epsilon_{k}$, where $\epsilon_{k}\overset{i.i.d.}{\sim}\mathcal{N}\,(0,\sigma_{\epsilon}^{2})$ for $1\leq k\leq K$. Following an EM algorithm developed in [11], the MAP estimate of $x_{k}$ given the observed data, $x_{k|K}$, is obtained. Then, the direct multitaper PSD of the estimated process $x_{k|K}$ is taken as the PSD estimate of $x_{k}$. For the multivariate non-stationary case, for each process $x_{k,j}$, $1\leq j\leq J$, we assume joint stationarity within non-overlapping consecutive windows of length $W$. The MAP estimate of each process $\mathbf{x}_{m,j}=[x_{(m-1)W+1,j},x_{(m-1)W+2,j},\cdots,x_{mW,j}]^{\top}$, for $1\leq m\leq M$ and $1\leq j\leq J$ are obtained, and an estimate of the ESD matrix is derived using the corresponding estimators. We refer to this estimator as the State-Space ESD (SS-ESD) estimator.

Peristimulus Time Histogram ESD Estimator

This estimator is derived by directly considering the ensemble mean of the spiking observations $\overline{n}_{k,j}$, referred to as the peristimulus time histogram (PSTH), as an estimate of the random signal $x_{k,j}$, for $1\leq k\leq K$ and $1\leq j\leq J$ [13]. With a similar joint stationarity assumption in windows of length $W$, the non-overlapping sliding window multitaper spectral estimate of the PSTH forms the PSTH-ESD estimator.

In order to benchmark our comparison, we consider the theoretical spectra of the AR processes derived using closed-form expressions (True ESD) as well as the the non-overlapping sliding window direct multitaper estimates of the processes $x_{k,j}$ that have been used to generate the spikes. We refer to the later benchmark as the Oracle ESD, as if an oracle could directly observe the latent processes and estimate their ESD.

V-A Case 1: Latent trivariate process observed through spiking

In order to simulate the three latent processes ($J=3$) with spectral couplings, we consider different linear combinations of a set of AR(6) processes, $\{\{y_{k}^{(i)}\}_{k=1}^{K},\,1\leq i\leq 6\}$, where $y_{k}^{(i)}$ is tuned around the frequency $f_{i}$, with $f_{1}=1.15$ Hz, $f_{2}=0.95$ Hz, $f_{3}=1.3$ Hz, $f_{4}=1.5$ Hz, $f_{5}=0.65$ Hz and $f_{6}=1.85$ Hz. All signals have been sampled at a sampling rate of $f_{s}=32$ Hz, for a total duration of $2000$ seconds ($K=64000$). The trivariate random process is defined as:

[TABLE]

where $x_{i,dc}$ are the DC components, $u_{k}$ is the unit step function, $\nu_{i,k}$ is a zero mean white Gaussian noise with unit variance, and $\sigma_{i}$ is a scaling standard deviation to set the SNR of all signals at $20$ dB, for $i=1,2$ and $3$.

In words, $x_{k,1}$ has been formed by combining three AR components, $y_{k}^{(1)}$, $y_{k}^{(2)}$ and $y_{k}^{(5)}$. The component related to $y_{k}^{(1)}$ has been amplitude modulated by a low frequency cosine signal at $f_{0}=0.0008$ Hz. The signal initially consists of only $y_{k}^{(1)}$ and $y_{k}^{(4)}$, and the third component $y_{k}^{(5)}$ is added to $x_{k,1}$ after $800$ seconds. The process $x_{k,2}$ is composed of four AR components, $y_{k}^{(2)}$, $y_{k}^{(4)}$, $y_{k}^{(5)}$ and $y_{k}^{(6)}$, with the contribution from $y_{k}^{(4)}$ having a lag of $6$ samples. The third process $x_{k,3}$ consists of the three AR components of $y_{k}^{(3)}$, $y_{k}^{(5)}$ and $y_{k}^{(6)}$. Although $y_{k}^{(6)}$ appears with a lag of $10$ samples initially, it becomes in phase with after $1000$ seconds.

We generate spike trains for $L=20$ neurons per each of the latent processes, using the logistic link model of Eq. (2). All DC components in the CIF have been set to $-5.5$, so that the average spiking rate of the ensemble corresponding to each CIF is $\approx 0.28$ spikes/second, consistent with the low spiking rate of experimentally recorded data. A $30$ second sample window of the process $x_{k,1}$ and the corresponding raster plot of the spiking observations is shown in Fig. 1.

We take the window length of stationary to be $100$ seconds ($W=3200$), resulting in a total number of $M=20$ windows. The parameter $\alpha$ has been fixed at $0.4$ to have optimal dependency across time windows and the prior parameter $\rho$ has been set to $0.2$. Note that these hyper-parameters can be further fine tuned using cross validation, in order to achieve a higher precision. Furthermore, we set $N=800$ in order to have a densely sampled spectral representation. Given that the spectral range is known to be $[0,2]$ Hz, we only use the first $N_{\max}=100$ frequency bins of the matrix $A$ for the sake of computational complexity. Finally, the time-bandwidth product of the multitaper framework is chosen as $2$ ($\xi=2$), and the first three dpss tapers are used.

Fig. 2 shows the estimation results corresponding to this case. The results have been formatted as a grid, with the columns representing (from left to right) the True EDS, Oracle ESD, PPMT-ESD, SS-ESD, and PSTH ESD estimates. The rows represent (from top to bottom) $(\bm{\Psi}_{m})_{1,1}(\omega)$, $(\bm{\Psi}_{m})_{2,2}(\omega)$, $(\bm{\Psi}_{m})_{3,3}(\omega)$, $(\bm{\Psi}_{m})_{1,2}(\omega)$, $(\bm{\Psi}_{m})_{2,3}(\omega)$ and $(\bm{\Psi}_{m})_{1,3}(\omega)$, where $(\bm{\Psi}_{m})_{i,j}(\omega)$ denotes the magnitude of the $(i,j)^{\sf th}$ block of the ESD matrix. Moreover, in order to have a closer inspection, the magnitude of the spectra corresponding to a window of $t=700$ s to $t=800$ s for $(\bm{\Psi}_{m})_{1,1}(\omega)$, $(\bm{\Psi}_{m})_{2,2}(\omega)$ and $(\bm{\Psi}_{m})_{1,2}(\omega)$ are depicted in Fig. 3.

It can be observed that the proposed PPMT-ESD estimator (Fig. 2(C)) results in much less background noise compared to all the others, while precisely capturing the dynamic evolution of the spectra and properly resolving the various frequency components. The latter is more evident from Fig. 3, where the PPMT-ESD (black trace) closely matches the true ESD (blue trace) on par with the Oracle ESD (red trace), while the SS-ESD (green trace) and PSTH-ESD (orange trace) show significant bias and variability. It is worth mentioning that the erroneous spectral peak near the DC component in Fig. 3 (black trace) is due to the estimation error of the DC component in the CIF model, in the low spiking regime of our setting. Nevertheless, this peak appears in the SS-ESD and PSTH-ESD estimates as well, but is mitigated by increasing the spiking rate.

Due to the time-domain smoothing model used in the SS-ESD estimator, the ESD rapidly decays with frequency (Fig. 3, green trace). As such, the SS-ESD estimate (Fig. 2(D)) heavily amplifies non-existing low frequency components that arise from the intrinsic noise in spiking observations, while suppressing the higher frequency components that exist in the true ESD (Fig. 2(A)). Similar, the estimate of the PSTH-ESD shown in Fig. 2(E) fails to capture most of the temporal and spectral features of the ESD, since it does not account for the binary nature of the observations. To quantify the performance of these estimators, we repeated this numerical experiment for a total of $50$ trials, generating independent realizations of the AR processes and spiking observations per trial, and computed the Mean Squared Error (MSE) with respect to the the True ESD (in dB scale). The average and variance of the MSE values are presented in Table I. All MSE computations have been normalized with respect to the total power of the True ESD (in dB scale). As expected, the Oracle ESD achieves the lowest MSE. Among the three methods that use the spiking observations, our proposed PPMT-ESD estimator achieves the lowest MSE, followed by the SS-ESD estimator with a significant gap. The PSTH-ESD estimator exhibits the poorest performance, in terms of both average MSE and variance, which is also visually evident in Fig. 3. In the spirit of easing reproducibility, a MATLAB implementation that regenerates the data, results and figures outlined in this section has been made publicly available on the open source repository GitHub [40].

V-B Case 2: Latent bivariate process with one directly observable component

While Case 1 was a natural choice for performance comparison, Case 2 is of particular interest in the joint analysis of neural spiking and continuous signals, such as the local field potential (LFP). The LFP signal corresponds to the electrical field potential measured at the cortical surface, and mesoscale dynamics of cortical activity.

In this case, we consider a bivariate random process, whose first component $x_{k,1}$ is observed through spiking activity $\{n_{k,1}^{(l)}\}_{k,l=1}^{K,L}$, while its second component $x_{k,2}$ is directly observable in i.i.d. zero-mean Gaussian noise, i.e., $\widetilde{x}_{k,2}:=x_{k,2}+\nu_{k}$, with $\nu_{k}\sim\mathcal{N}\,(0,\sigma_{\nu}^{2})$. Using the same notations as before, the two processes used in this simulation are,

[TABLE]

The process $x_{k,2}$ considered in this scenario is exactly the same as that described in Case 1. Process $x_{k,1}$ is a slightly modified version of the process $x_{k,1}$ in Case 1: the first two components of $x_{k,1}$ are the same as described before; the amplitude modulated component related to $y_{k}^{(1)}$ and the fixed frequency AR component $y_{k}^{(4)}$. However, instead of the third component $y_{k}^{(5)}$, we include a frequency modulated component, $y_{k}^{(7)}$, which is an AR process tuned around the frequency $f_{7}$. The frequency $f_{7}$ changes by decrements of $0.06$ Hz every $200$ seconds, starting at $0.9$ Hz at $t=0$ s.

The ESD matrix of this bivariate process can be estimated by Algorithm 2 with a slight modification in the Forward filtering step (step 2) of Algorithm 1: Given that the second process is directly observable, the distribution $f(\mathcal{D}_{1}^{m}|\{\mathbf{V}\}_{1}^{m},\widehat{\bm{\theta}}^{(r)})$ needs to be modified, and accordingly, the log-posterior in Eq. (13) changes to:

[TABLE]

The parameters used in forming the estimates are the same as those chosen for Case 1. Fig. 4 shows the estimation results corresponding to this case. The results have been similarly formatted as a grid, with the columns representing (from left to right) the True EDS, Oracle ESD, PPMT-ESD, SS-ESD, and PSTH ESD estimates. The rows represent (from top to bottom) $(\bm{\Psi}_{m})_{1,1}(\omega)$, $(\bm{\Psi}_{m})_{2,2}(\omega)$, and $(\bm{\Psi}_{m})_{1,2}(\omega)$. Note that in this case, we take the SS-ESD and PSTH-ESD estimates of $(\bm{\Psi}_{m})_{2,2}(\omega)$ to be the same as its Oracle ESD estimate, given that these methods are based on estimating the process $x_{k,2}$ in time domain (which is directly observable here).

Similar to Case 1, the proposed PPMT-ESD estimator (Fig. 4(C)) captures the dynamics of the spectra $(\bm{\Psi}_{m})_{1,1}(\omega)$ and $(\bm{\Psi}_{m})_{1,2}(\omega)$ accurately, closely matching the True ESD (Fig. 4(A)). As before, the SS-ESD estimator (Fig. 4(D)) is not able to capture the ESD dynamics, especially at high frequencies. Though some frequency components at certain time windows are detected by the PSTH-ESD estimates (Fig. 4(E)), most of the frequency content is concealed by background noise.

VI Application to experimentally-recorded data from general anesthesia

Finally, we apply our proposed PPMT-ESD estimator to multi-unit recordings from a human subject under Propofol-induced general anesthesia (data from [12, 14]). The data set includes the spiking activity of $41$ neurons as well as the LFP recorded from a patient undergoing intra-cranial monitoring for surgical treatment of epilepsy using a multichannel micro-electrode array implanted in temporal cortex [12]. Recordings were conducted during the administration of Propofol for induction of anesthesia. The experimental protocol is extensively explained in [12]. Since the original recordings have been over-sampled, we down-sample both the LFP and spike recordings to the same sampling frequency of $25$ Hz [14].

From the spike recordings, we select $25$ neurons ($L=25$) with the highest spiking rates for analysis. The average spiking rate of this subpopulation of neurons is $0.1732$ spikes/second. Fig. 5(a) shows the LFP signal for a total duration of $800$ s considered in the analysis. A zoomed-in view of the LFP and the raster plot of the neuronal ensemble corresponding to a $40$ s window from $t=480$ s to $t=520$ s is shown in Fig. 5.

Similar to Case 2 in Section V-B, we model the neuronal ensemble and the LFP by a bivariate process, where the LFP is directly observable, but the latent process driving the neuronal activity is observed through spiking. We assume the process to be stationary within windows of length $40$ s ($W=1000$), resulting in $M=20$ non-overlapping windows. Given that the main spectral content is known to be in the range of [math]–$2$ Hz [12, 14], we estimate the ESDs up to $2.5$ Hz ($N_{\max}=100$). Further, we choose $\xi=2$, and use the first three tapers in the analysis. We set $\alpha=0.85$ and $\rho=0.02$.

Fig. 6 shows the results of our analysis. The first column illustrates the ESD estimates, where the dashed vertical line marks the induction of anesthetic. The second column shows snapshots of the normalized magnitude spectrum corresponding to a $40$ s window from $t=480$ s to $t=520$ s. The first three rows (from top to bottom) show the estimates of $(\bm{\Psi}_{m})_{1,1}(\omega)$, i.e., the ESD of the latent process driving spiking activity, for PSTH-ESD, SS-PSD, and the proposed PPMT-ESD, respectively. The fourth row represents the Oracle ESD estimate of the LFP signal ($(\bm{\Psi}_{m})_{1,1}(\omega)$), and the fifth row depicts ($(\bm{\Psi}_{m})_{1,2}(\omega)$), the cross-spectra between the LFP and the latent process, using the proposed PPMT-ESD method.

Similar to our simulation studies, we observe that the PSTH-ESD estimate (Fig. 6, first row) captures a significant number of spurious frequency components and harmonics, which are known to be absent during Propofol general anesthesia [12], masking the relevant spectral content. Also, as before, the SS-ESD estimate amplifies the low-frequency components (Fig. 6, second row). Even though the dominant spectral content of the LFP is around $0.5$ Hz and $1.1$ Hz, as observed in the right panel of the fourth row in Fig. 6, neither of these peaks are present in the SS-ESD estimate (Fig. 6, second row, right panel). In addition, the rapid change in the LFP spectrogram (fourth row, left panel) marked by the dashed line, is not captured by the SS-ESD estimate (second row, left panel).

The proposed PPMT-ESD estimate (third row) closely resembles the spectrum of the LFP (fourth row), evident in both the spectrogram and the spectral snapshot. The dominant frequency components of the PMTM-ESD around $0.5$ Hz and $1.1$ Hz match those in the LFP spectrum. The temporal variations of the spectrum are also consistent with that of the LFP, as visible in the left panels. The estimated cross-spectrum (fifth row) further corroborates this observation, by capturing the spectral coupling between the LFP and the latent process driving the spiking activity. It is indeed hypothesized that the latent process driving the spiking activity is the LFP signal itself [12]. As such, our analysis corroborates this hypothesis by providing an accurate estimate of the ESDs of both the LFP and the latent process driving the neuronal spiking.

VII Concluding Remarks

Brain oscillations are known to play a significant role in modulating the finescale dynamics of neuronal spiking. These oscillations are also known to be non-stationary, as they reflect the changes in the internal states of behavioral conditions. Understanding how the latent processes that drive spiking activity are related to brain oscillations is a key problem in computational neuroscience. On one hand, the theory of point processes has been successfully utilized in recent years to capture the dynamics of neuronal spiking data. On the other hand, spectral estimation of non-stationary continuous signals has been widely studied. Existing methods either treat the spiking data as continuous, and apply non-stationary spectral estimation techniques off-the-shelf, or try to first estimate these latent processes in time-domain, followed by forming spectrograms. Both these approaches are known to suffer from high bias and variability. A unified methodology for inferring spectrotemporal representations of the multivariate latent non-stationary processes that drive neuronal spiking is lacking. In this work, we proposed such a methodology for estimating the ESD matrix of a multivariate non-stationary latent process directly from binary spiking observations. To this end, we integrated techniques from state-space modeling, multitaper analysis and point processes.

We established theoretical bounds on the bias and variance performance of the proposed estimator, and compared its performance with the aforementioned existing techniques through application to simulated and experimentally-recorded neural data. Our simulation studies confirmed our theoretical analysis and revealed the favorable performance of our proposed method over existing approaches. Our application to real data corroborated the hypothesis on the role of the local field potential in regulating spiking activity under general anesthesia, by providing a clear picture of the underlying spectral couplings. While we have developed our methodology in the context of neuronal data analysis, it can be applied to a wide range of discrete observations that are modulated by underlying oscillations, such as heart-beat data [41]. Our methodology can also be extended to infer non-stationary network-level properties such as the frequency domain Granger-Geweke causality [42, 43].

VIII Acknowledgment

We would like to thank Emery N. Brown and Patrick L. Purdon for sharing the data from [12].

Appendix A Proof of Theorem 1

Proof.

Let $S(\omega)$ be the PSD of the process $\{x_{k}\}_{k=1}^{K}$. Then,

[TABLE]

where $(a)$ follows from the triangle inequality. Further,

[TABLE]

where $(b)$ follows from the Cauchy-Schwarz inequality. Thus, the desired bounds on the bias and variance can be established through bounding the first and second moments of $(\widehat{S}^{\sf mt}(\omega)\,-\,S^{\sf mt}(\omega))$. The first moment can be bounded by

[TABLE]

[TABLE]

where $(c)$ and $(d)$ follow from the triangle inequality and $(e)$ follows by bounding the complex sinusoid. The main technical difficulty in further development of the bounds is due to the fact that $\operatorname{logit}(z)$ does not have a Taylor series expansion for $z\in(0,1)$. We thus have to find other algebraically useful bounds. To this end, we need to establish the following technical lemma.

Lemma 1**.**

Consider the event $A=\{\overline{n}_{k}\;|\;\overline{n}_{k}\neq 0,\,\overline{n}_{k}\neq 1,\,1\leq k\leq K\}$. The following inequality holds true for all $\overline{n}_{k}\in A$:

[TABLE]

where

[TABLE]

Proof of Lemma 1.

First, consider the case $\lambda_{k}\leq 0.5$. We bound the function $\varepsilon(\overline{n}_{k})$ in a piece-wise fashion as follows. Note that $\operatorname{logit}(\overline{n}_{k})$ is convex for $\overline{n}_{k}\geq 0.5$ and concave for $\overline{n}_{k}\leq 0.5$. Thus, it immediately follows that for $\overline{n}_{k}\leq\lambda_{k}$, $\varepsilon(\overline{n}_{k})$ is convex and hence,

[TABLE]

Furthermore, for $\lambda_{k}\leq\overline{n}_{k}\leq 0.5$, $\varepsilon(\overline{n}_{k})$ is concave, and hence is bounded by the tangent at $\lambda_{k}$ as follows.

[TABLE]

Finally, for the case of $\overline{n}_{k}\geq 0.5$, consider the line

[TABLE]

From the convexity of $\varepsilon(\overline{n}_{k})$, $\ell(\overline{n}_{k})$ upper bounds $\varepsilon(\overline{n}_{k})$ for $\overline{n}_{k}\geq 0.5$, since $\ell(0.5)\geq\varepsilon(0.5)$ for $\lambda_{k}\leq 0.5$. Combining the piece-wise bounds in Eqs. (26), (27) and (28), we conclude the claim in Lemma 1 for $\lambda_{k}\leq 0.5$. Due to the symmetry of $\varepsilon(\overline{n}_{k})$, through a similar argument, the bound can be established for $\lambda_{k}>0.5$, which concludes the proof. ∎

Given that $|x_{k}|\;\leq\,B$ and assuming that $L$ is large enough so that $L\geq 2(1+\exp(B))$, we can further simplify the bound of Lemma 1. We have:

[TABLE]

Thus, for sufficiently large $L$, we conclude that,

[TABLE]

Now, consider the expectations in the bounds of Eq. (25). Using iterated conditioning,

[TABLE]

where $(f)$ follows from triangle and Jensen’s inequalities. In order further simplify these bounds, we invoke the result of Lemma 1. First, note that $\operatorname{logit}(\overline{n}_{k})$ is unbounded in the complement of event $A$. Provided $|x_{k}|\,\leq\,B$, $\mathbb{P}(\overline{n}_{k}\neq 0)$ and $\mathbb{P}(\overline{n}_{k}\neq 1)$ can be lower bounded by $1-\exp(-L\log(1+\exp(-B)))$, which implies that

[TABLE]

Therefore, for sufficiently large $L$, we see that $\mathbb{P}(A)$ is exponentially close to $1$. Thus, hereafter we condition the expectations on the highly probable event $A$. From Eq. (29), we get

[TABLE]

Note that the random variable $n_{k}=L\overline{n}_{k}$ is the sum of $L$ independent Bernoulli random variables given $x_{k}$. Thus, given $x_{k}$, $n_{k}\sim\operatorname{Binomial}(L,\lambda_{k})$. Accordingly,

[TABLE]

where $(g)$ follows from the Jensen’s inequality and $(h)$ follows from substituting expression for the variance of a binomial random variable. Further, note that $\lambda_{k}(1-\lambda_{k})\leq 1/4$, for $\lambda_{k}\in[0,1]$ and $\mathbb{P}(A)=1-(\lambda_{k}^{L}+(1-\lambda_{k})^{L})\geq 1/2$, if $(1+\exp(B))^{L}>2(1+\exp(BL))$, which is satisfied for large enough $L$. Thus, combining the bounds in Eqs. (31) and (29), we get,

[TABLE]

By a similar argument we can show that,

[TABLE]

Thus, the expectation in Eq. (30) is bounded as:

[TABLE]

Following a similar argument, one can show for $n\neq m$,

[TABLE]

Finally, using the bounds of Eqs. (33) and (34) and noting that $\sum_{k=0}^{K-1}(\nu_{k}^{(p)})^{2}=1$, for all $0\leq p\leq P$, we can upper bound the expectation in Eq. (25) as,

[TABLE]

where

[TABLE]

This concludes the proof of the bound on bias. Following similar bounding techniques, the second moment in Eq. (24) can be bounded by:

[TABLE]

where

[TABLE]

This concludes the proof of Theorem 1. ∎

Appendix B Proof of Corollaries 1 and 2

Proof of Corollary 1.

Proof of Corollary 1 follows the proof of Theorem 1 closely, with the natural extension to the multivariate case. Following the proof of Theorem 1, the constants $g^{\prime}_{1}$ and $g^{\prime}_{2}$ in this case are given by:

[TABLE]

[TABLE]

∎

Proof of Corollary 2.

As for Corollary 2, we work under the technical assumptions of Theorems 1 and 2 in [25]. Following [25], we assume that in Eq. (10), $\mathbf{Q}_{m}=\mathbf{Q}$ for all $m$ in this proof, and that the EM algorithm finds estimates of $\mathbf{Q}$ and $\alpha$ close to their true value (for large enough $K$). Then, under Assumptions (1) and (2), we identify the effective observation $\widetilde{\mathbf{y}}_{m}^{(p)}$ corresponding to the $p^{\sf th}$ taper at window $m$ in [25] by the concatenation of \nu_{w}^{(p)}\operatorname{logit}\big{(}\bar{n}_{(m-1)W+w,j}\big{)} for $w=1,2,\cdots,W$ and $j=1,2,\cdots,J$ in a vector of length $WJ$. We also assume, without loss of generality that $W=uN$, for some integer $u$. Then, we denote by $\bm{\Sigma}_{\infty}$ the steady state covariance of the backward smoother, and $\bm{\Lambda}:=\alpha\bm{\Sigma}_{\infty}(\alpha^{2}\bm{\Sigma}_{\infty}+\mathbf{Q})^{-1}$ and $\bm{\Gamma}=(\alpha^{2}\bm{\Sigma}_{\infty}+\mathbf{Q})[\mathbf{I}-uW((\alpha^{2}\bm{\Sigma}_{\infty}+\mathbf{Q})^{-1}+uW\mathbf{I})^{-1}]$, as in [25]. Note that these matrices are $NJ\times NJ$ in our case. Under the same assumptions [25], we consider them to be diagonal with the $i^{\sf th}$ diagonal element being $\gamma_{i}$ and $\eta_{i}$ respectively. Then, following the proof of Theorem 1 and those of Theorems 1 and 2 in [25], it can be shown that the statement of the corollary holds with the constants:

[TABLE]

[TABLE]

and

[TABLE]

∎