Estimating summary statistics in the spike-train space

Abstract

Estimating sample averages and sample variability is important in analyzing neural spike trains data in computational neuroscience. Current approaches have focused on advancing the use of parametric or semiparametric probability models of the underlying stochastic process, where the probabilistic distribution is characterized at each time point with basic statistics such as mean and variance. To directly capture and analyze the average and variability in the observation space of the spike trains, we focus on a data-driven approach where statistics are defined and computed in a function space in which the spike trains are viewed as individual points. Based on the definition of a “Euclidean” metric, a recent paper introduced the notion of the mean of a set of spike trains and developed an efficient algorithm to compute it under some restrictive conditions. Here we extend this study by: (1) developing a novel algorithm for mean computation that is quite general, and (2) introducing a notion of covariance of a set of spike trains. Specifically, we estimate the covariance matrix using the geometry of the warping functions that map the mean spike train to each of the spike trains in the dataset. Results from simulations as well as a neural recording in primate motor cortex indicate that the proposed mean and covariance successfully capture the observed variability in spike trains. In addition, a “Gaussian-type” probability model (defined using the estimated mean and covariance) reasonably characterizes the distribution of the spike trains and achieves a desirable performance in the classification of the spike trains.

Appendices

Appendix A: Existence of the optimal warping function

We at first state one theoretical fact on the optimal time warping on the metric d _p , which is actually the same as Lemma 1 in Wu and Srivastava (2011). We restate it here to emphasize its importance.

Lemma 1

For p ≥ 1 and any positive diffeomorphism γ from [a, c] to [b, d] (i.e. $\gamma(a) = b, \gamma(c) = d, 0 < \dot \gamma(t) < \infty$ ), the time warping has the following optimal cost:

$$\inf_{\gamma \in \Gamma} \int_a^c \left|1 - \dot \gamma(t)^{1/p} \right|^p dt = \left|(c-a)^{1/p} - (d-b)^{1/p}\right|^p. $$

The infimum is reached if the warping is linear from [a, c] to [b, d] (i.e. $\dot \gamma(t) = \frac{d-b}{c-a}$ ). This condition is also necessary for any p > 1.

Based on this result, we have the following theorem on the existence of the optimal warping function between any two spike trains.

Theorem 1

Assume $f(t) = \sum_{i=1}^M \delta(t-t_i) \in {\cal S}_M$ and $g(t) = \sum_{j=1}^N \delta(t-s_j) \in {\cal S}_N$ are two spike trains in [0, T]. Then there exists a warping function $\gamma^\ast \in \Gamma_{f,g}$ , such that

$$ \begin{array}{rll} d_p[\lambda] (f,g) &=& \left(X([f], [g \circ \gamma^\ast])\vphantom{\int_{0}^T}\right.\\ &&\ \ \ \left.+ \lambda \int_{0}^T \left|1 - \dot \gamma^\ast(t)^{1/p}\right|^p dt \right)^{1/p}\ , \end{array} $$

(14)

where 1 ≤ p < ∞ and 0 < λ < ∞.

Proof

For any warping function γ ∈ Γ, we have two sets of spike times, [f] and [g ∘ γ]. Let $\{a_l\}_{l=1}^L$ be the intersection of [f] and [g ∘ γ], where L is the cardinality of {a _l }. This implies X([f],[g ∘ γ]) = M + N − 2L. Denoting a ₀ = 0 and a _L + 1 = T, we define a warping function $\tilde \gamma \in \Gamma_{f,g}$ as follows:

$$ \begin{array}{rll} \tilde \gamma(t) &=& \frac{\gamma(a_{l}) -\gamma(a_{l-1})}{a_{l}-a_{l-1}} (t - a_{l-1}) \\&&+ \gamma(a_{l-1}), \mbox { if } a_{l-1} \le t \le a_l, l = 1, \cdots, L+1. \end{array} $$

For illustration, we show the warping functions γ and $\tilde \gamma$ in Fig. 12, where $\tilde \gamma$ is basically a piecewise linear version of γ by linearly connecting all grid points (a _l , γ(a _l )), l = 0, ⋯ , L + 1. We note that $\tilde \gamma$ remains the time matchings between {a _l } and {γ(a _l )}. This implies that

$$X([f],[g\circ \tilde \gamma]) \le M + N - 2L = X([f],[g\circ \gamma]). $$

As $\tilde \gamma$ is linear from [a _l − 1, a _l ] to [γ(a _l − 1), γ(a _l )], using Lemma 1, we have

$$ \begin{array}{rll} &&\int_{a_{l-1}}^{a_l} \left|1 - {\dot {\tilde \gamma}(t)}^{1/p}\right|^p dt \\ && \quad\le \int_{a_{l-1}}^{a_l} \left|1 - {\dot { \gamma}(t)}^{1/p}\right|^p dt, \ \ \ l = 1, \cdots, L+1. \end{array} $$

Fig. 12

Any warping function γ ∈ Γ (green dashed line) and a function $\tilde \gamma \in \Gamma_{f,g}$ (blue solid line) obtained by linearly interpolating γ at the intersection of [f] and [g ∘ γ]

Therefore,

$$ \begin{array}{rll} && X([f],[g\circ \tilde \gamma]) +\lambda \int_0^T \left|1 - {\dot {\tilde \gamma}(t)}^{1/p}\right|^p dt \\ &&\quad \le X([f],[g\circ \gamma]) + \lambda \int_0^T \left|1 - {\dot { \gamma}(t)}^{1/p}\right|^p dt. \end{array} $$

This indicates $\tilde \gamma \in \Gamma_{f,g}$ provides a better time warping to measure the distance between f and g than γ. Note that the set Γ _f,g is finite. This implies that,

$$ \begin{array}{rll} &&\inf\limits_{\gamma \in \Gamma} \left ( X([f],[g\circ \gamma]) +\lambda \int_0^T \left|1 - {\dot { \gamma}(t)}^{1/p}\right|^p dt \right )^{1/p} \\ &&\quad = \min\limits_{\gamma \in \Gamma_{f,g}} \left (X([f],[g\circ \gamma]) + \lambda \int_0^T \left|1 - {\dot { \gamma}(t)}^{1/p}\right|^p dt \right )^{1/p} . \end{array} $$

Therefore, there exists $\gamma^\ast \in \Gamma_{f, g} \subset \Gamma$, such that

$$ \begin{array}{rll} d_p[\lambda] (f,g) &=& \left(X([f], [g \circ \gamma^\ast])\vphantom{\int_{0}^T}\right.\\ &&\ \ \ \ \left.+ \lambda \int_{0}^T \left|1 - {\dot \gamma}^\ast(t)^{1/p}\right|^p dt \right)^{1/p}\ . \end{array} $$

□

Appendix B: Convergence of the MCP-algorithm

Theorem 2

Denote the estimated mean in the j th iteration of the MCP-algorithm as S ^(j) . Then the sum of squared distances $\sum_{i=1}^N d_2(S_i, S^{(j)})^2$ decreases iteratively. That is,

$$\sum\limits_{i=1}^N d_2\left(S_i, S^{(j+1)}\right)^2 \le \sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\right)^2.$$

Proof

In the jth iteration, we at first perform the Matching Step and found the optimal warping γ _i from S _i to S ^(j), i = 1, ⋯ , N. That is,

$$ d_2\left(S_i, S^{(j)}\right)^2 = X\left(\left[S_i \circ \gamma_i\right], \left[S^{(j)}\right]\right) + \lambda \left|\left|1-\sqrt{\dot \gamma_i} \right|\right|^2. $$

Then, we perform the Centering Step to compute the extrinsic mean of warping functions $\{\gamma_i\}_{i=1}^N$. As shown in Section 2.2, the mean warping function, $\bar \gamma$, of $\{\gamma_i\}_{i=1}^N$ is also piecewise linear and its SRVF has the closed form given by Eq. (6). By the definition of the extrinsic mean, we have

$$ \begin{array}{rll} \sum\limits_{i=1}^N \left|\left|1 - \sqrt{\gamma_i}\right|\right|^2 &\ge& \sum\limits_{i=1}^N \left|\left| \sqrt{\dot {\bar \gamma}} - \sqrt{\dot \gamma_i}\right|\right|^2 \\[6pt] &=& \sum\limits_{i=1}^N \left|\left|1 - \sqrt{ {\gamma_i \dot \circ {\bar \gamma}^{-1}}}\right|\right|^2 \end{array} $$

We then apply ${\bar \gamma}^{-1}$ to both S _i  ∘ γ _i and S ^(j). In this operation, the spikes in both trains are moved with the same warping function. Using the definition of the d ₂ metric, we have

$$ \begin{array}{rll} &&\sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\circ {\bar \gamma}^{-1}\right)^2 \\ && \quad\le \sum\limits_{i=1}^NX\left(\left[S_i \circ \left(\gamma_i \circ {\bar \gamma}^{-1}\right)\right], \left[S^{(j)} \circ {\bar \gamma}^{-1}\right]\right) \\ &&\qquad +\, \lambda \left|\left|1-\sqrt{ \gamma_i \dot \circ {\bar \gamma}^{-1}} \right|\right|^2. \end{array} $$

One can easily verify that

$$ X\left(\left[S_i \circ \gamma_i\right], \left[S^{(j)}\right]\right) = X\left(\left[S_i \circ \gamma_i \circ {\bar \gamma}^{-1}\right], \left[S^{(j)} \circ {\bar \gamma}^{-1}\right]\right). $$

Therefore,

$$ \begin{array}{rll} \sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\circ {\bar \gamma}^{-1}\right)^2&\le& \sum\limits_{i=1}^N X\left(\left[S_i \circ \gamma_i\right], \left[S^{(j)}\right]\right) \\ && +\, \lambda \left|\left|1-\sqrt{\dot \gamma_i} \right|\right|^2 \\ &=& \sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\right)^2. \end{array} $$

Next, we perform the Pruning Step and update the number of spikes in the mean. $[\widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}]$ is the set of spikes in $[S^{(j)} \circ {\bar \gamma^{-1}}]$ that appear more than N/2 times in {[$S_i\circ\gamma_i \circ \bar\gamma^{-1}$]}. Base on the definition of the d ₂ metric, we have

$$ \begin{array}{rll} &&\sum\limits_{i=1}^N d_2\left(S_i, \widetilde {S^{(j)}\circ {\bar \gamma}^{-1}}\right)^2 \\ &&\quad \le \sum\limits_{i=1}^NX\left(\left[S_i \circ \left(\gamma_i \circ {\bar \gamma}^{-1}\right)\right], \left[\widetilde {S^{(j)} \circ {\bar \gamma}^{-1}}\right]\right)\\ &&\qquad +\, \lambda \left|\left|1-\sqrt{ \gamma_i \dot \circ {\bar \gamma}^{-1}} \right|\right|^2. \end{array} $$

Using the basic rule of the Exclusive OR, it is easy to find that

$$ \begin{array}{rll} &&\sum\limits_{i=1}^N X\left(\left[S_i \circ \left(\gamma_i \circ {\bar \gamma}^{-1}\right)\right], \left[\widetilde {S^{(j)} \circ {\bar \gamma}^{-1}}\right]\right) \\ &&\quad \le \sum\limits_{i=1}^N X\left(\left[S_i \circ \left(\gamma_i \circ {\bar \gamma}^{-1}\right)\right], \left[S^{(j)} \circ {\bar \gamma}^{-1}\right]\right). \end{array} $$

Therefore,

$$ \begin{array}{rll} &&\sum\limits_{i=1}^N d_2\left(S_i, \widetilde {S^{(j)}\circ {\bar \gamma}^{-1}}\right)^2\\ &&\quad\le \sum\limits_{i=1}^NX\left(\left[S_i \circ \left(\gamma_i \circ {\bar \gamma}^{-1}\right)\right], \left[ {S^{(j)} \circ {\bar \gamma}^{-1}}\right]\right) \\ &&\qquad +\, \lambda \left|\left|1-\sqrt{ \gamma_i \dot \circ {\bar \gamma}^{-1}} \right|\right|^2 \le \sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\right)^2 \end{array} $$

Finally, we test whether one spike in $[\widetilde {S^{(j)}\circ {\bar \gamma}^{-1}}]$ that appears the minimal number of times in {[$S_i\circ\gamma_i \circ \bar\gamma^{-1}$]} can be removed.

• (i) If the spike can be removed, the remaining set of spikes is denoted as $[\widetilde {S^{(j)}\circ {\bar \gamma}^{-1}}^{-}]$. This is under the condition that $\sum_{i=1}^N d_2(S_i, \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}^{-})^2 \le \sum_{i=1}^N d_2(S_i, \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}})^2$. In this case, we update the mean spike train $S^{(j+1)} = \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}^{-}$, and the number of spikes in the mean n = |S ^(j + 1)|. In this case, we have

  $$ \begin{array}{rll} \sum\limits_{i=1}^N d_2\left(S_i, S^{(j+1)}\right)^2 &=& \sum\limits_{i=1}^N d_2\left(S_i, \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}^{-}\right)^2\\ &\le& \sum\limits_{i=1}^N d_2\left(S_i, \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}\right)^2 \\ &\le& \sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\right)^2 \end{array} $$

• (ii) If the spike cannot he removed, we update the mean spike train $S^{(j+1)} = \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}$, and the number of spikes in the mean n = |S ^(j + 1)|. In this case, we also have

  $$ \begin{array}{rll} \displaystyle\sum\limits_{i=1}^N d_2\left(S_i, S^{(j+1)}\right)^2 &=& \sum\limits_{i=1}^N d_2\left(S_i, \widetilde {S^{(j)} \circ {\bar \gamma^{-1}}}\right)^2\\ &\le& \sum\limits_{i=1}^N d_2\left(S_i, S^{(j)}\right)^2. \end{array} $$

In summary, we have shown that the sum of squared distances $\sum_{i=1}^N d_2(S_i, S^{(j)})^2$ decreases iteratively.□