Feature extraction from spike trains with Bayesian binning: ‘Latency is where the signal starts’

Abstract

The peristimulus time histogram (PSTH) and its more continuous cousin, the spike density function (SDF) are staples in the analytic toolkit of neurophysiologists. The former is usually obtained by binning spike trains, whereas the standard method for the latter is smoothing with a Gaussian kernel. Selection of a bin width or a kernel size is often done in an relatively arbitrary fashion, even though there have been recent attempts to remedy this situation (DiMatteo et al., Biometrika 88(4):1055–1071, 2001; Shimazaki and Shinomoto 2007a, Neural Comput 19(6):1503–1527, 2007b, c; Cunningham et al. 2008). We develop an exact Bayesian, generative model approach to estimating PSTHs. Advantages of our scheme include automatic complexity control and error bars on its predictions. We show how to perform feature extraction on spike trains in a principled way, exemplified through latency and firing rate posterior distribution evaluations on repeated and single trial data. We also demonstrate using both simulated and real neuronal data that our approach provides a more accurate estimates of the PSTH and the latency than current competing methods. We employ the posterior distributions for an information theoretic analysis of the neural code comprised of latency and firing rate of neurons in high-level visual area STSa. A software implementation of our method is available at the machine learning open source software repository (www.mloss.org, project ‘binsdfc’).

Appendices

Appendix A: Computing the evidence with dynamic programming

The evidence, or marginal likelihood of a model with M bins is given by (see Eq. (9)):

$$\label{app_evidence} \begin{array}{rcl} P(\{\vec{z}^i\}|M)&=&\sum\limits_{k_{M-1}=M-1}^{T-2} \sum\limits_{k_{M-2}=M-2}^{k_{M-1}-1} \ldots \\ && \ldots \sum\limits_{k_0=0}^{k_1-1} P(\{\vec{z}^i\}|\{k_m\},M)P(\{k_m\}|M) \end{array} $$

(31)

where the summation boundaries are chosen such that the bins are non-overlapping and contiguous and

$$\label{app_pgivenk} \begin{array}{l} P(\{\vec{z}^i\}|\{k_m\},M)\\ {\kern20pt} =\int_0^1 d\{f_m\} P(\{\vec{z}^i\} | \{f_m\},\{k_m\},M) p(\{f_m\}|M). \end{array} $$

(32)

Recall that the probability of a (multi)set of spiketrains $\{\vec{z}^i\}=\{z_1,\ldots,z_N\}$, assuming independent generation, is given by Eq. (3):

$$\label{app_likelihood_set} \begin{array}{rcl} P(\{\vec{z}^i\}|\{f_m\},\{k_m\},M)&=&\prod\limits_{i=1}^N \prod\limits_{m=0}^M f_m^{s(\vec{z}^i,m)}(1-f_m)^{g(\vec{z}^i,m)} \\ &=& \prod\limits_{m=0}^M f_m^{s(\{\vec{z}^i\},m)}(1-f_m)^{g(\{\vec{z}^i\},m)} \\ \end{array} $$

(33)

where $s(\{\vec{z}^i\},m)=\sum_{i=1}^N s(\vec{z}^i,m)$ is the number of spikes in all spiketrains in bin m and $g(\{\vec{z}^i\},m)=\sum_{i=1}^N g(\vec{z}^i,m)$ is the number of all non-spikes, or gaps. The prior of the firing rates (Eq. (5)) is

$$ \label{app_fprior} p(\{f_m\}|M) = \prod_{m=0}^M \mbox{B}(f_m;\sigma_m,\gamma_m). $$

(34)

The integrals in Eq. (32) can be evaluated by virtue of Eqs. (33 and 34):

$$ \label{pzk} P(\{\vec{z}^i\}|\{k_m\},M)\!=\!\prod_{m=0}^M \! \frac{\mbox{B}(s(\{\vec{z}^i\},m)\!+\!\sigma_m,g(\{\vec{z}^i\},m)+\gamma_m)}{\mbox{B}(\sigma_m,\gamma_m)} $$

(35)

where $\mbox{B}(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$ is Euler’s Beta function (Davis 1972). Equation (35) is a product with one factor per bin, and each factor depends only on spike/gap counts and prior parameters in that bin. To compute Eq. (31), we can therefore use an approach very similar to that described in (Endres and Földiák 2005; Endres 2006) in the context of density estimation and in (Hutter 2006, 2007) for Bayesian function approximation: define the function

$$\label{intervalevidences_def} \begin{array}{rcl} \mbox{getIEC}({k_s},{k_e},{m})&:=&\mbox{B}(s(\{\vec{z}^i\},k_s,k_e)\\ &&+\sigma_m,g(\{\vec{z}^i\},k_s,k_e)+\gamma_m) \end{array} $$

(36)

where $s(\{\vec{z}^i\},k_s,k_e)$ is the number of spikes and $g(\{\vec{z}^i\},k_s,k_e)$ is the number of gaps in $\{\vec{z}^i\}$ between the start interval k _s and the end interval k _e (both included). Furthermore, collect all contributions to Eq. (31) that do not depend on the data (i.e. $\{\vec{z}^i\}$) and store them in the array $\mbox{pr}[M]$:

$$ \label{priors_def} \mbox{pr}[M]:=\frac{\prod_{m=0}^M \frac{1}{\mbox{B}(\sigma_m,\gamma_m)}}{\left(\begin{array}{c}T-1\\M\end{array}\right)}. $$

(37)

Substituting Eq. (35) into Eq. (31) and using the definitions (36) and (37), we obtain

$$\label{evidence_redef} \begin{array}{rcl} P(\{\vec{z}^i\}|M) &\propto& \sum\limits_{k_{M-1}=M-1}^{T-2}\ldots \\ &&\ldots \sum\limits_{k_0=0}^{k_1-1} \prod\limits_{m=1}^M \mbox{getIEC}({k_{m-1}+1},{k_m},{m}) \\ && \times \mbox{getIEC}({0},{k_0},{0}) \end{array} $$

(38)

with k _M  = T − 1 and the constant of proportionality being $\mbox{pr}[M]$. Since the factors on the r.h.s. depend only on two consecutive bin boundaries each, it is possible to apply dynamic programming (Bertsekas 2000): rewrite the r.h.s. by ‘pushing’ the sums as far to the right as possible:

$$\label{pushed_sums} \begin{array}{rcl} &&{\kern-7pt} P(\{\vec{z}^i\}|M)\\ &&{\kern12pt} \propto \sum\limits_{k_{M-1}=M-1}^{T-2} \mbox{getIEC}({k_{M-1}\!+\!1},{T\!-\!1},{M}) \\ &&{\kern20pt} \times \sum\limits_{k_{M-2}=M-2}^{k_{M-1}-1} \mbox{getIEC}({k_{M-2}+1},{k_{M-1}},{M-1}) \\ &&{\kern20pt} \times \ldots \sum\limits_{k_0=0}^{k_1-1} \mbox{getIEC}({k_0+1},{k_1},{1})\mbox{getIEC}({0},{k_0},{0}).\\ \end{array} $$

(39)

Evaluating the sum over k ₀ requires O(T) operations (assuming that T ≫ M, which is likely to be the case in real-world applications). As the summands depend also on k ₁, we need to repeat this evaluation O(T) times, i.e. summing out k ₀ for all possible values of k ₁ requires O(T ²) operations. This procedure is then repeated for the remaining M − 1 sums, yielding a total computational effort of O(MT ²). Thus, initialise the array $\mbox{subE}_{0}[k]:=\mbox{getIEC}({0},{k},{0})$, and iterate for all m = 1,...,M:

$$ \mbox{subE}_{m}[k]:=\sum_{r=m-1}^{k-1} \mbox{getIEC}({r+1},{k},{m}) \mbox{subE}_{m-1}[r], $$

(40)

A close look at Eq. (39) reveals that while we sum over k _M − 1, we need $\mbox{subE}_{M-1}[k]$ for k = M − 1;...;T − 2 to compute the evidence of a model with its latest boundary at T − 1. We can, however, compute $\mbox{subE}_{M-1}[T-1]$ with little extra effort, which is, up to a factor $\mbox{pr}[M-1]$, equal to $P(\{\vec{z}^i\}|M-1)$, i.e. the evidence for a model with M − 1 bin boundaries. Moreover, having computed $\mbox{subE}_{k}[m]$, we do not need $\mbox{subE}_{m-1}[k-1]$ anymore. Hence, the array $\mbox{subE}_{k}[m-1]$ can be reused to store $\mbox{subE}_{k}[m]$, if overwritten in reverse order. Table 4 shows this algorithm in pseudo-code ($\mbox{E}[m]$ contains the evidence of a model with m bin boundaries inside [t _min,t _max] after termination).

Table 4 Computing the evidences of models with up to M bin boundaries

Appendix B: Computing the posterior distribution of the latency

We compute the joint probability of the latency L = t and the observed spiketrains $\{\vec{z}^i\}$ given the number of bins and the signal separation level S via

$$ \begin{array}{rcl} P(L=t,\{\vec{z}^i\}|M,S)&=&\sum\limits_{k_{M-1}=M-1}^{T-2} \sum\limits_{k_{M-2}=M-2}^{k_{M-1}-1} \ldots \\ && \ldots \sum\limits_{k_0=0}^{k_1-1} P(L=t,\{\vec{z}^i\},\{k_m\}|M,S) \end{array} $$

(41)

where

$$ \begin{array}{rcl} P(L&=&t,\{\vec{z}^i\},\{k_m\}|M,S) \\ &=&\int_0^1 d\{f_m\} P(L=t|\{k_m\},\{f_m\},M,S) \\ &&\times\,p(\{\vec{z}^i\},\{f_m\},\{k_m\}|M). \end{array} $$

(42)

Note that P(L = t|{k _m },{f _m },M,S) is the r.h.s of Eq. (21) and $p(\{\vec{z}^i\},\{f_m\},\{k_m\}|M)$ is the numerator of the r.h.s. of Eq. (8). As a consequence of Eq. (21), the only nonzero contributions to the average are models which have a (lower) bin boundary at t. Assume t was at the lower bound of bin j, i.e. at t = k _j − 1 + 1 (the {k _m } are inclusive upper bin boundaries, as defined above). Carrying out the integrals over the {f _m } yields:

$$\label{latency_intbounds} \begin{array}{rcl} P(L&=&t,\{\vec{z}^i\},\{k_m\}|M,S) \\ &=& \int_0^1 d\{f_m\} P(L=t|\{k_m\},\{f_m\},M,S) \\ &&\times\, p(\{\vec{z}^i\},\{f_m\},\{k_m\}|M) \\ &=& \int_0^S df_0 \ldots \int_0^S df_{j-1} \int_S^1 df_{j} \int_0^1 df_{j+1} \ldots \\ &&\ldots\,\int_0^1 df_M p(\{\vec{z}^i\},\{f_m\},\{k_m\}|M) \end{array} $$

(43)

The integration boundaries in the last line of Eq. (43) are a consequence of our latency definition: all bins m < j will contribute to the integral only as long as f _m  < S, hence the upper bound of their integrals is at S. Bin j contributes only if f _j  ≥ S, thus the lower bound of the integral over f _j is S. The bins m > j are not affected by the latency probability (Eq. (21)) whence their integrals still run from 0 to 1. By virtue of Eqs. (3) and (5), we obtain

$$\label{latency_pzk} \begin{array}{rcl} P(L&=&t,\{\vec{z}^i\},\{k_m\}|M,S) \\ &=& \prod\limits_{m=0}^{j-1} B_S(s(\{\vec{z}^i\},m)+\sigma_m,g(\{\vec{z}^i\},m)+\gamma_m) \\ &&\times\, \tilde{B}_S(s(\{\vec{z}^i\},j)+\sigma_m,g(\{\vec{z}^i\},j)+\gamma_m)) \\ &&\times\, \prod\limits_{m=j+1}^M B(s(\{\vec{z}^i\},m)+\sigma_m,g(\{\vec{z}^i\},m)+\gamma_m) \\ &&\times\, \prod\limits_{m=0}^M \frac{1}{B(\sigma_m,\gamma_m)} P(\{k_m\}|M) \end{array} $$

(44)

where $B_S(a,b)=\int_0^S t^{a-1}(1-t)^{b-1}dt$ is the incomplete Beta function (Davis 1972) and $\tilde{B}_S(a,b)=\int_S^1 t^{a-1}(1-t)^{b-1}dt=B(a,b)-B_S(a,b)$ is its complement. Note that up to the factor P({k _m }|M), Eq. (44) is basically Eq. (35) with some of the Beta functions replaced by incomplete Beta functions. Hence, the remaining summations over the {k _m } can be carried out by using

$$\label{intervallatencies_def} \begin{array}{rcl} && \mbox{getIEC}_{L}({k_s},{k_e},{m}):= \\ &&{\kern15pt} =\left\{\begin{array}{ll} B_S(s(\{\vec{z}^i\},k_s,k_e) & \\[6pt] {\kern6pt}+\sigma_m,g(\{\vec{z}^i\},k_s,k_e)+\gamma_m)&\,\,\,\mbox{if } k_e<t\\[6pt] \tilde{B}_S(s(\{\vec{z}^i\},k_s,k_e) & \\[6pt] {\kern6pt}+\sigma_m,g(\{\vec{z}^i\},k_s,k_e)+\gamma_m))&\,\,\,\mbox{if } k_s=t\\[6pt] B(s(\{\vec{z}^i\},k_s,k_e) & \\[6pt] {\kern6pt} +\sigma_m,g(\{\vec{z}^i\},k_s,k_e)+\gamma_m)&\,\,\,\mbox{if } k_s>t \\[6pt] 0 & \,\,\,\mbox{otherwise}\end{array}\right. \end{array} $$

(45)

instead of $\mbox{getIEC}({k_s},{k_e},{m})$ (Eq. (36)) in the evidence computation algorithm. This procedure yields the desired $P(L=t,\{\vec{z}^i\}|M,S)$. The above derivation is an instance of the general framework for computing expectations of functions of bin boundaries and firing probabilities described in Endres and Földiák (2005).

Appendix C: Computing the posterior density of the firing rate

We compute the joint probability density of the firing rate $f(t) = \tilde{f}$, the latency L = t and the observed spiketrains $\{\vec{z}^i\}$ given the number of bins and the signal separation level S via

$$ \begin{array}{rcl} &&{\kern-7pt} P(f(t) = \tilde{f},L=t,\{\vec{z}^i\}|M,S) = \\ &&{\kern12pt} =\sum\limits_{k_{M-1}=M-1}^{T-2} \sum\limits_{k_{M-2}=M-2}^{k_{M-1}-1} \ldots \\ &&{\kern11pt} \ldots \sum\limits_{k_0=0}^{k_1-1} P(f(t) = \tilde{f},L=t,\{\vec{z}^i\},\{k_m\}|M,S) \end{array} $$

(46)

where

$$ \begin{array}{rcl} &&{\kern-7pt} p(f(t) = \tilde{f},L=t,\{\vec{z}^i\},\{k_m\}|M,S) \\ &&{\kern12pt} =\int_0^1 d\{f_m\} P(f(t) = \tilde{f},L=t|\{k_m\},\{f_m\},M,S) \\ &&{\kern20pt} \times\, p(\{\vec{z}^i\},\{f_m\},\{k_m\}|M) \end{array} $$

(47)

Note that $P(f(t) = \tilde{f},L=t|\{k_m\},\{f_m\},M,S)$ is the r.h.s of Eq. (24) and $p(\{\vec{z}^i\},\{f_m\},\{k_m\}|M)$ is the numerator of the r.h.s. of Eq. (8). As a consequence of Eq. (24), the only nonzero contributions to the average are models which have a (lower) bin boundary at t and f(t) ≥ S. Assume t was at the lower bound of bin j, i.e. at t = k _j − 1 + 1 (the {k _m } are inclusive upper bin boundaries, as defined above). Integrating out the {f _m } yields

$$\label{latfire_pzk} \begin{array}{rcl} &&{\kern-7pt} p(f(t)= \tilde{f},L=t,\{\vec{z}^i\},\{k_m\}|M,S) \\ &&{\kern12pt} =\prod\limits_{m=0}^{j-1} B_S(s(\{\vec{z}^i\},m)+\sigma_m,g(\{\vec{z}^i\},m)+\gamma_m) \\ &&{\kern20pt} \times \tilde{f}^{s(\{\vec{z}^i\},j)+\sigma_m-1}\;(1-\tilde{f})^{g(\{\vec{z}^i\},j)+\gamma_m)-1} \\ &&{\kern20pt} \times \prod\limits_{m=j+1}^M B(s(\{\vec{z}^i\},m)+\sigma_m,g(\{\vec{z}^i\},m)+\gamma_m) \\ &&{\kern20pt} \times \prod\limits_{m=0}^M \frac{1}{B(\sigma_m,\gamma_m)} P(\{k_m\}|M) \end{array} $$

(48)

where the second line is a result of Eq. (3) and (5) multiplied with the Dirac delta function in Eq. (23). Hence, the remaining summations over the {k _m } can be carried out by using

$$\label{intervallatfire_def} \begin{array}{rcl} &&{\kern-7pt} \mbox{getIEC}_{f,L}({k_s},{k_e},{m}):= \\ &&{\kern12pt} \left\{\begin{array}{ll} B_S(s(\{\vec{z}^i\},k_s,k_e) & \\[6pt] {\kern6pt}+\sigma_m,g(\{\vec{z}^i\},k_s,k_e)+\gamma_m)&\,\,\,\mbox{if } k_e<t \\[6pt] \tilde{f}^{s(\{\vec{z}^i\},k_s,k_e)+\sigma_m-1} & \\[6pt] {\kern6pt}\times (1-\tilde{f})^{g(\{\vec{z}^i\},k_s,k_e)+\gamma_m -1}&\,\,\,\mbox{if } k_s=t \\[6pt] B(s(\{\vec{z}^i\},k_s,k_e) & \\[6pt] {\kern6pt}+\sigma_m,g(\{\vec{z}^i\},k_s,k_e)+\gamma_m)&\,\,\,\mbox{if } k_s>t \\[6pt] 0 & \,\,\,\mbox{otherwise}\end{array}\right. \end{array} $$

(49)

instead of $\mbox{getIEC}({k_s},{k_e},{m})$ (Eq. (36)) in the evidence computation algorithm. Thus we obtain $P(f(t) = \tilde{f},L=t,\{\vec{z}^i\}|M,S)$.