An information-geometric framework for statistical inferences in the neural spike train space

Abstract

Statistical inferences are essentially important in analyzing neural spike trains in computational neuroscience. Current approaches have followed a general inference paradigm where a parametric probability model is often used to characterize the temporal evolution of the underlying stochastic processes. To directly capture the overall variability and distribution in the space of the spike trains, we focus on a data-driven approach where statistics are defined and computed in the function space in which spike trains are viewed as individual points. To this end, we at first develop a parametrized family of metrics that takes into account different warpings in the time domain and generalizes several currently used spike train distances. These new metrics are essentially penalized L ^p norms, involving appropriate functions of spike trains, with penalties associated with time-warping. The notions of means and variances of spike trains are then defined based on the new metrics when p = 2 (corresponding to the “Euclidean distance”). Using some restrictive conditions, we present an efficient recursive algorithm, termed Matching-Minimization algorithm, to compute the sample mean of a set of spike trains with arbitrary numbers of spikes. The proposed metrics as well as the mean spike trains are demonstrated using simulations as well as an experimental recording from the motor cortex. It is found that all these methods achieve desirable performance and the results support the success of this novel framework.

Appendices

Appendix A: Metric D _p [σ, λ]

Here we need to show the distance D _p [σ, λ] satisfies the following three conditions in the definition of a metric space, where p ≥ 1, σ > 0, and λ > 0.

(a)  Non-negativity: D(f,g) ≥ 0 for any f, g $\in \mathcal{S}^{(\sigma)}$. The equality holds if and only if f = g.

Proof

The non-negativity is apparent.□

(b)  Symmetry: D(f,g) = D(g,f) for any f, g $\in \mathcal{S}^{(\sigma)}$.

Proof

For any γ(t) ∈ Γ, let s = γ(t). Then t = γ ^− 1(s), and s = γ(γ ^− 1(s)), $1 = \dot \gamma (\gamma^{-1}(s)) \dot \gamma^{-1}(s)$. Therefore, $\int |f(t)^{1/p}$ − $[g(\gamma(t))\dot \gamma(t)]^{1/p}|^p$ + λ|1 − $\dot \gamma(t)^{1/p}|^p dt$ = $\int |[f(\gamma^{-1}(s))\dot \gamma^{-1}(s)]^{1/p}$ − g(s)^1/p|^p + $\lambda |\dot \gamma^{-1}(s)^{1/p}$ − 1|^p ds. Taking infimum over γ(t) and γ ^− 1(s) on both sides, the symmetry holds.□

(c)  Triangle-inequality: D(f,g) ≤ D(f,h) + D(h,g) for any f, g, h $\in \mathcal{S}^{(\sigma)}$.

Proof

The triangle-inequality will hold if for any γ(t), γ ₁(t) ∈ Γ we have

$$ \begin{array}{rll} &&\left[\int \!\left|f(t)^{1/p} - \left[g(\gamma(t))\dot \gamma(t)\right]^{1/p}\right|^p\! + \lambda \left|1 - \dot \gamma(t)^{1/p}\right|^p dt\right]^{1/p} \!\!\!\!\!\!\!\!\!\!\!\!\nonumber \\ &&{\kern3pt}\le \left[\int \left|f(t)^{1/p} - \left[h(\gamma_1(t))\dot \gamma_1(t)\right]^{1/p}\right|^p\right.\nonumber\\ &&\qquad+\left. \lambda \left|1 - \dot \gamma_1(t)^{1/p}\right|^p dt\vphantom{\int}\right]^{1/p} \nonumber \\ &&{\kern10pt} +\, \left[\int \left|h(s)^{1/p} - \left[g(\gamma \circ \gamma_1^{-1}(s)) \gamma \dot \circ \gamma_1^{-1}(s)\right]^{1/p}\right|^p \right.\nonumber\\ &&{\kern4pt}\qquad +\left. \lambda \left|1 - \gamma \dot \circ \gamma_1^{-1}(s)^{1/p}\right|^p ds\vphantom{\int}\right]^{1/p} \label{eq:tiq} \end{array} $$

(14)

where ∘ denotes function composition in Γ (and therefore $\gamma \circ \gamma_1^{-1}$ is still a warping in the group Γ). Let s = γ ₁(t), the second term on the right hand side of Eq. (14) is equal to

$$ \begin{array}{rll} &\left[\int \left|\left[h(\gamma_1(t))\dot \gamma_1(t)\right]^{1/p} - \left[g(\gamma(t))\dot \gamma(t)\right]^{1/p}\right|^p \right.\\ &{\kern6pt} +\left. \lambda \left|\dot \gamma_1(t)^{1/p} - \dot \gamma(t)^{1/p}\right|^p dt\vphantom{\int}\right]^{1/p} \end{array} $$

Therefore, Eq. (14) will hold if we can show that for any a, b, c, d, e, and f in the L ^p space

$$ \begin{array}{rll} & \left[\int |a - b|^p + |d - e|^p dt\right]^{1/p} \notag\\ &\le \left[\int |a - c|^p + |d - f|^p dt\right]^{1/p} \notag\\ &{\kern6pt} + \left[\int |c - b|^p + |f - e|^p dt\right]^{1/p} \label{eq:tiq2} \end{array} $$

(15)

Indeed, Eq. (15) is a direct result from the triangle inequalities in the L ^p space and the L ^p norm in the Euclidean $\Re^2$ vector space.□

Appendix B: Metric d _p [λ]

(B.1)  At first, we show the distance d _p [λ] also satisfies the following three conditions in the definition of a metric space, where p ≥ 1 and λ > 0.

• (a) Non-negativity: d(f,g) ≥ 0 for any f, g $\in \mathcal{S}$. The equality holds if and only if f = g.

Proof

The non-negativity is apparent.□

• (b) Symmetry: d(f, g)  =  d(g,f) for any f, g $\in\! \mathcal{S}$.

Proof

Assume $f \in \cal{S}_M$ and $g \in \cal{S}_N$. For any γ(t) ∈ Γ, let s = γ(t). Then t = γ ^− 1(s), and s = γ(γ ^− 1(s)), $1 = \dot \gamma (\gamma^{-1}(s)) \dot \gamma^{-1}(s)$. Therefore, $M \!+\! N \!-\! 2\sum_{i=1}^M \sum_{j=1}^N {\bf 1}_{t_i = \gamma(s_j)} \!+ \lambda \int |1 - \dot \gamma(t)^{1/p}|^p dt =$ $N \!+\! M \!-\! 2\sum_{j=1}^N \sum_{i=1}^M {\bf 1}_{s_j = \gamma^{-1}(t_i)} \!+ \lambda \int |\dot \gamma^{-1}(s)^{1/p} -$ 1|^p ds. Taking infimum over γ(t) and γ ^− 1(s) on both sides, the symmetry holds.□

• (c) Triangle-inequality: d(f,g) ≤ d(f,h) + d(h, g) for any f, g, h $\in \mathcal{S}$.

Proof

Assume $f \! = \! \sum_{l=i}^L \delta(t \! - \! r_i)$, $g \! = \! \sum_{j=1}^M \delta(t -$ s _j ), and h = $\sum_{k=1}^N \delta(t$ − t _k ). The triangle-inequality will hold if for any γ(t), γ ₁(t) ∈ Γ we have

$$ \begin{array}{rll} && \left[L + M - 2\sum\limits_{i=1}^L \sum\limits_{j=1}^M {\bf 1}_{r_i = \gamma(s_j)} \right.\nonumber\\ &&+\left. \lambda \int \left|1 - \dot \gamma(t)^{1/p}\right|^p dt\vphantom{\sum\limits_{i=1}^L}\right]^{1/p} \nonumber\\ &&\le \left[L + N - 2\sum\limits_{i=1}^L \sum\limits_{k=1}^N {\bf 1}_{r_i = \gamma_1(t_k)} \right.\nonumber\\ &&\,\,\,\quad+ \left.\lambda \int \left|1 - \dot \gamma_1(t)^{1/p}\right|^p dt\vphantom{\sum\limits_{k=1}^N}\right]^{1/p} \label{eq:tiqq} \nonumber\\ &&\,\,\,\,+ \left[N + M - 2\sum\limits_{k=1}^N \sum\limits_{j=1}^M {\bf 1}_{t_k = \gamma_1^{-1} \circ \gamma(s_j)} \right.\nonumber\\[-5pt] &&\,\,\,\,\,\,\,\quad+\left. \lambda \int |1 - \dot {\gamma_1^{-1} \circ \gamma}(t)^{1/p}|^p dt\vphantom{\sum\limits_{j=1}^M}\right]^{1/p} \end{array} $$

(16)

where ∘ denotes function composition in Γ. Let s = γ(t) in the first and third integrations and let s = γ ₁(t) in the the second integration, Eq. (16) can be written as

$$ \begin{array}{rll} &&\left[L + M - 2\sum\limits_{i=1}^L \sum\limits_{j=1}^M {\bf 1}_{r_i = \gamma(s_j)} \right.\\ &&{\kern3pt}+ \left.\lambda \int |1 - \dot {\gamma^{-1}}(s)^{1/p}|^p ds\vphantom{\sum\limits_{i=1}^L}\right]^{1/p} \\ &&{\kern3pt}\le \left[L + N - 2\sum_{i=1}^L \sum_{k=1}^N {\bf 1}_{r_i = \gamma_1(t_k)} \right.\\ &&\,\,\qquad + \left.\lambda \int |1 - \dot {\gamma_1^{-1}}(s)^{1/p}|^p ds\vphantom{\sum\limits_{i=1}^L}\right]^{1/p}\\ &&{\kern10pt}+ \left[N + M - 2\sum\limits_{k=1}^N \sum\limits_{j=1}^M {\bf 1}_{\gamma_1(t_k) = \gamma(s_j)} \right.\\ &&\qquad\quad+\left. \lambda \int |\dot {\gamma^{-1}}(s)^{1/p} - \dot {\gamma_1^{-1}}(s)^{1/p}|^p ds\vphantom{\sum\limits_{i=1}^L}\right]^{1/p} \end{array} $$

It is straightforward to verify that

$$ \begin{array}{rll} N + \sum\limits_{i=1}^L \sum\limits_{j=1}^M {\bf 1}_{r_i = \gamma(s_j)} &\ge& \sum\limits_{i=1}^L \sum\limits_{k=1}^N {\bf 1}_{r_i = \gamma_1(t_k)} \\ &&+\, \sum\limits_{k=1}^N \sum\limits_{j=1}^M {\bf 1}_{\gamma_1(t_k) = \gamma(s_j)}. \end{array} $$

This implies that

$$ \begin{array}{rll} && L + M - 2\sum\limits_{i=1}^L \sum\limits_{j=1}^M {\bf 1}_{r_i = \gamma(s_j)}\\ &&{\kern3pt}\le L + N - 2\sum\limits_{i=1}^L \sum\limits_{k=1}^N {\bf 1}_{r_i = \gamma_1(t_k)}\\ &&{\kern9pt} +\, N + M - 2\sum\limits_{k=1}^N \sum\limits_{j=1}^M {\bf 1}_{\gamma_1(t_k) = \gamma(s_j)}. \end{array} $$

Therefore, Eq. (16) will hold if we can show that for any non-negative values x, y, z with x ≤ y + z and a, b, c in the L ^p space

$$ \begin{array}{rll} \left[\!x +\! \int |a - b|^p dt\!\right]^{1/p} \!\le{}& \!\left[y + \!\int |a - c|^p dt\right]^{1/p} \nonumber\\ &+\! \left[z + \!\int\! |c - b|^p dt\right]^{1/p} \label{eq:tiqq2} \end{array} $$

(17)

Indeed, Eq. (17) is also a direct result from the triangle inequalities in the L ^p space and the L ^p norm in the Euclidean $\Re^2$ vector space.□

(B.2)  Then, we show that when σ → 0, the limiting form of D _p [σ, λ] is d _p [λ].

Proof

Assume $f(t) \!=\! \sum_{i=1}^M K_{\sigma}(t\!-\!t_i) \!\in\! {\cal S}_M^{(\sigma)}$ and $g(t) = \sum_{j=1}^N K_{\sigma}(t-s_j) \in {\cal S}_N^{(\sigma)}$ are two smoothed spike trains in [0, T]. When σ is sufficiently small, spikes in f and g will either fully match (two spike kernels are entirely overlapped) or do not match at all (two spike kernels are well separated). Moreover, when two spike kernels K _σ (t − t _i ) in f and K _σ (t − s _j ) in g are fully matched, there must be no time warping in this match. That is, γ(s _j ) = t _i and $\dot \gamma (t) = 1$ when s _j  − σ < t < s _j  + σ.

This indicates that

$$ \begin{array}{rll} &\left|{f(t)}^{1/p} - [g(\gamma(t))\dot \gamma(t)]^{1/p}\right|^p \notag\\ &\,\,\,= \left \{ \begin{array}{@{}l@{\,\,}l} 0 & \mbox{ if fully matched} \\ f(t) + g(\gamma(t))\dot \gamma(t) & \mbox{ if not matched} \end{array} \right . \end{array} $$

(18)

Therefore, when σ is sufficiently small, $\int_0^T |$ ${f(t)}^{1/p} \!-\! [g(\gamma(t))\dot \gamma(t)]^{1/p}|^p dt = M + N - 2\sum_{i=1}^M$ $\sum_{j=1}^N {\bf 1}_{t_i = \gamma(s_j)}$. This is actually the total number of spike kernels that are not overlapped after warping. As the penalty term in d _p [λ] is exactly the same as that in D _p [σ, λ], d _p [λ] is the limiting form of D _p [σ, λ] when σ→0.□

Appendix C: Limiting form for upper bounds

Here we show that the upper bound is exact in some limiting form for both D ₂ and D ₁ distances. Our proof is based on an example plot in Fig. 11, which shows an infinite sequence of piecewise linear warpings from [0, T] to [0, T]. Computing the integral of the piecewise linear warping, we have

$$ \begin{array}{rll} \int_0^T \sqrt{\dot \gamma(t)} dt &= \int_0^{\frac{n-1}{n}T} \sqrt{\frac{1}{n-1}} dt + \int_{\frac{n-1}{n}T}^T \sqrt{n-1} dt \\ &= \frac{2T}{n}\sqrt{n-1} \rightarrow 0 (n \rightarrow \infty) \end{array} $$

and

$$ \begin{array}{rll} \int_0^T \! |1 \! - \! \dot \gamma(t)| dt \! &=& \! \int_0^{\frac{n-1}{n}T} \! \left( \! 1 \! - \! \frac{1}{n \! - \! 1} \! \right) dt \! + \! \int_{\frac{n -1}{n}T}^T (n \! - \! 1 \! - \! 1) dt \\ &=& \frac{2(n-2)}{n}T \rightarrow 2T (n \rightarrow \infty). \end{array} $$

This proves that the upper bound, 2 λT, is exact for both D ₂ and D ₁.

Fig. 11

A piecewise linear warping from [0, T] to [0, T]

Appendix D: Optimal time warping

Here we state some theoretical facts on the optimal time warping.

Lemma 1

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

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

When p > 1, the equality holds if and only if the warping is linear from [a, c] to [b, d]. When p = 1, such linear property is sufficient.

Proof

In domain [a, c], for any function f(t) ∈ L ^p, we denote its L ^p norm as $||f(t)||_p = (\int_a^c |f(t)|^{p}| dt)^{1/p}$. Then,

$$ \begin{array}{rll} \int_a^c \left|1 - {\dot \gamma(t)}^{1/p}\right|^p dt &=& \left\| 1 - {\dot \gamma(t)}^{1/p} \right\|_p^p \\ &\ge& \left| \| 1 \|_p - \left\| {\dot \gamma(t)}^{1/p} \right\|_p\right | ^p \\ &=& \left|(c-a)^{1/p} - (d-b)^{1/p}\right|^p. \end{array} $$

Using theory in the L ^p spaces, when p > 1, the equality holds if and only if $\dot \gamma(t)$ is a positive constant, i.e. the warping is linear from [a, c] to [b, d]. When p = 1, the equality holds if the sign of $1 - \dot \gamma(t)$ does not change in [a, c] (always nonnegative or always nonpositive), for example, if the warping is linear from [a, c] to [b, d].□

Theorem 1

Let spike trains $f, g \in \mathcal{S}_M$ in [0, T] and the corresponding ISIs are Δs ₁, ⋯ , Δs _M + 1 and Δt ₁, ⋯ , Δt _M + 1 , respectively. If λ < 1/(2^p − 1 T), then

$$ d_p[\lambda](f,g) = \left (\lambda \sum\limits_{k=1}^{M+1}|(\Delta s_k)^{1/p} - (\Delta t_k)^{1/p}|^p \right )^{1/p} \label{eq:lm2dis} $$

(19)

Proof

In domain [0, T], for any function h(t) ∈ L ^p, we denote its L ^p norm as $||f(t)||_p = (\int_0^T |f(t)|^{p} dt)^{1/p}$. Then,

$$ \begin{array}{rll} \int_0^T \left|1 - {\dot \gamma(t)}^{1/p}\right|^p dt &=& \left\| 1 - {\dot \gamma(t)}^{1/p} \right\|_p^p \\ &\le& \left| \| 1 \|_p + \left\| {\dot \gamma(t)}^{1/p} \right\|_p\right | ^p = 2^pT. \end{array} $$

This indicates that the total penalty is bounded by λ2^p T. If one spike in one train does not match with any spike in the other train, the cost in the matching term should be at least 2. Therefore, when λ2^p T < 2, or λ < 1/(2^p − 1 T), all M spikes in f and all M spikes in g must be one-to-one matched for the optimal distance.

When all spikes are matched, the distance between two spike trains only depends on the optimal warping cost for each pair of ISIs. The distance in Eq. (19) is a direct result from Lemma 1.□

Theorem 2

Let spike trains $f \in \mathcal{S}_M$ and $g \in \mathcal{S}_N$ in [0, T]. If λ < 1/(2^p − 1 T), then

$$ \begin{array}{rll} &d_p[\lambda](f,g) \notag\\ &\,\,\,= \left (|M-N|+\lambda \cdot (\mbox{optimal warping cost in the}\right.\notag\\ &\quad\qquad\qquad\qquad\qquad\left. \mbox{associated ISIs}) \right )^{1/p}. \end{array} $$

Proof

Without loss of generality, we assume M ≥ N. The matching cost is minimal with value M − N when all N spikes in g are matched with N (of M) spikes in f. In this case, we have N + 1 pairs of ISIs associated the N spikes in f and N spikes in g. As shown in Theorem 1, the total penalty is always bounded by λ2^p T. Therefore, when λ < 1/(2^p − 1 T), the N spikes in g must match N (of M) spikes in f and the optimal total cost which results in the following distance:

$$ \begin{array}{rll} &d_p[\lambda](f,g) \\ &\,\,\,= \left (M-N+\lambda \cdot (\mbox{optimal warping cost in the $N+1$}\right.\\ &\quad\qquad\qquad\qquad\quad\left. \mbox{pairs of ISIs}) \right )^{1/p}. \end{array} $$

Moreover, using Theorem 1, if the N + 1 pairs of ISIs are $\{\Delta s_k\}_{k=1}^{N+1}$ and $\{\Delta t_k\}_{k=1}^{N+1}$, then the optimal time warping cost is: $\sum_{k=1}^{N+1} |{\Delta s_k}^{1/p} - {\Delta t_k}^{1/p}|^p.$ Therefore, the distance between f and g can be further represented in the following form:

$$ d_p[\lambda](f,g) = \left(\!(M - N) + \lambda \sum\limits_{k=1}^{J+1} |{\Delta s_k}^{1/p} - {\Delta t_k}^{1/p}|^p\!\right )^{1/p}. $$

□

Appendix E: Convergence of the MM-algorithm

Denote the estimated mean in the ith iteration of the MM-algorithm as S ⁽ⁱ⁾. Now we show that the sum of squared distance $\sum_{k=1}^K d_2^2(S_k, S^{(i)})$ decreases iteratively. That is,

$$ \sum_{k=1}^K d_2^2(S_k, S^{(i)}) \le \sum_{k=1}^K d_2^2(S_k, S^{(i-1)}). $$

As 0 is a natural lower bound, the iteration will always converge.

Proof

In the ith iteration, we have the initial S ^{(i − 1)} and let d _{k,i − 1} denote the distance between S _k and S ^{(i − 1)} based on current spike train matching subset. The number of spikes in the mean, n, is known beforehand. Hence the squared distance on the matching part is always |n _k  − n|, which is invariant with respect to any subset of spikes being used in the matching process. Therefore, we only need to focus on how the optimal matching cost changes from S ^{(i − 1)} to S ⁽ⁱ⁾.

In case 1, we compute the distance, r _k , between S _k and S ^{(i − 1)} and find n + 1 ISIs in S _k that have the optimal time warping with the n + 1 ISIs in S ^{(i − 1)}. As the ISIs are optimally selected, r _k  ≤ d _{k,i − 1}. In case 2, The descrease in distance also holds. However, as n _k  < n, the number of ISIs is only n _k  + 1, which could not be used to compute the mean using the closed-form formula in Eq. (11). In the MM-algorithm, we propose to linearly interpolate n − n _k imaginary spikes in S _k based on the n _k pairs of spikes between S _k and S ^{(i − 1)}. Here we show this interpolation does not change the warping distance between S _k and S ^{(i − 1)}.

Indeed, as the warping distance is represented in the sum of squared form (see the proof in Theorem 2), we only need to show the invariance in each squared term. Let Δa be one ISI in S _k and Δb be the associated ISI in S ^{(i − 1)}. Their squared distance is

$$ \left(\sqrt {\Delta a} - \sqrt {\Delta b}\right)^2. $$

Assume there are h spikes in the ISI in S ^{(i − 1)}, which partition Δb into h + 1 ISIs, denoted by (u ₁, ⋯ , u _h + 1) with $\sum_{j=1}^{h+1} u_j = \Delta b$. Our algorithm inserts h spikes in the ISI in S _k using linear interpolation. Therefore the (h + 1) ISIs in S _k have lengths

$$ \frac{\Delta a}{\Delta b} (u_1, \cdots, u_{h+1}). $$

Then the new warping distance based on these ISIs are

$$ \begin{array}{rll} \sum\limits_{j=1}^{h+1}\left(\sqrt {u_j} - \sqrt {\frac{\Delta a}{\Delta b}u_j}\right)^2 &=& \sum\limits_{j=1}^{h+1} u_j \left(1 - \sqrt {\frac{\Delta a}{\Delta b}}\right)^2 \\ &=& \Delta b \left(1 - \sqrt {\frac{\Delta a}{\Delta b}}\right)^2 \\ &=& \left(\sqrt {\Delta a} - \sqrt {\Delta b}\right)^2 \end{array} $$

In summary, over step 2, the distance between S _k and S ^{(i − 1)} has decreased and we have found a set of n + 1 ISIs in S _k that correspond to the n + 1 ISIs in S ^{(i − 1)}. In step 3, we update the ISIs in S ^{(i − 1)} to S ⁽ⁱ⁾. This update is based on the closed-form formula in Eq. (11) which minimizes the sum of squared distance. Therefore,

$$ \sum r_k^2 \ge \sum d_{k,i}^2. $$

As r _k  ≤ d _{k,i − 1}, k = 1, ⋯ , K, we have

$$ \sum d_{k,i-1}^2 \ge \sum d_{k,i}^2. $$

□