Mathematical methods

We define a compact formalism for multiunit spike-train metrics by opportunely characterising the space of multiunit feature vectors as a tensor product. Previous results from Houghton and Kreuz (2012, On the efficient calculation of van Rossum distances. Network: Computation in Neural Systems, 2012, 23, 48-58) on a clever formula for multiunit Van Rossum metrics are then re-derived within this framework, also fixing some errors in the original calculations.

A compact formalism for kernel-based multiunit spike train metrics

Consider a network with \(C\) cells. Let

\[\mathcal{U}= \left\{ \boldsymbol{u}^1, \boldsymbol{u}^2, \ldots, \boldsymbol{u}^C \right \}\]

be an observation of network activity, where

\[\boldsymbol{u}^i = \left\{ u_1^i, u_2^i, \ldots, u_{N_{\boldsymbol{u}^i}}^i \right \}\]

is the (ordered) set of times of the spikes emitted by cell \(i\). Let \(\mathcal{V}= \left\{\boldsymbol{v}^1, \boldsymbol{v}^2, \ldots, \boldsymbol{v}^C\right\}\) be another observation, different in general from \(\mathcal{U}\).

To compute a kernel based multiunit distance between \(\mathcal{U}\) and \(\mathcal{V}\), we map them to the tensor product space \(\mathcal{S} \doteq \mathbb{R}^C\bigotimes L_2(\mathbb{R}\rightarrow\mathbb{R})\) by defining

\[|\mathcal{U}\rangle = \sum_{i=1}^C |i\rangle \otimes |\boldsymbol{u}^i\rangle\]

where we consider \(\mathbb{R}^C\) and \(L_2(\mathbb{R}\rightarrow\mathbb{R})\) to be equipped with the usual euclidean distances, consequently inducing an euclidean metric structure on \(\mathcal{S}\) too.

Conceptually, the set of vectors \(\left\{ |i\rangle \right\}_{i=1}^C \subset \mathbb{R}^C\) represents the different cells, while each \(|\boldsymbol{u}^i\rangle \in L_2(\mathbb{R}\rightarrow\mathbb{R})\) represents the convolution of a spike train of cell \(i\) with a real-valued feature function \(\phi: \mathbb{R}\rightarrow\mathbb{R}\),

\[\langle t|\boldsymbol{u}\rangle = \sum_{n=1}^{N_{\boldsymbol{u}}}\phi(t-u_n)\]

In practice, we will never use the feature functions directly, but we will be only interested in the inner products of the \(|i\rangle\) and \(|\boldsymbol{u}\rangle\) vectors. We call \(c_{ij}\doteq\langle i|j \rangle=\langle i|j \rangle_{\mathbb{R}^C}=c_{ji}\) the multiunit mixing coefficient for cells \(i\) and \(j\), and \(\langle \boldsymbol{u}|\boldsymbol{v}\rangle=\langle \boldsymbol{u}|\boldsymbol{v}\rangle_{L_2}\) the single-unit inner product,

\[\begin{split}\label{eq:singleunit_intprod} \begin{split} \langle \boldsymbol{u}|\boldsymbol{v}\rangle & = \langle \left\{ u_1, u_2, \ldots, u_{N}\right\}|\left\{ v_1, v_2, \ldots, v_{M}\right\} \rangle = \\ &= \int\textrm{d }\!t \langle \boldsymbol{u}|t \rangle\langle t|\boldsymbol{v}\rangle = \int\textrm{d }\!t \sum_{n=1}^N\sum_{m=1}^M\phi\left(t-u_n\right)\phi\left(t-v_m\right)\\ &\doteq \sum_{n=1}^N\sum_{m=1}^M \mathcal{K}(u_n,v_m) \end{split}\end{split}\]

where \(\mathcal{K}(t_1,t_2)\doteq\int\textrm{d }\!t \left[\phi\left(t-t_1\right)\phi\left(t-t_2\right)\right]\) is the single-unit metric kernel, and where we have used the fact that the feature function \(\phi\) is real-valued. It follows immediately from the definition above that \(\langle \boldsymbol{u}|\boldsymbol{v}\rangle=\langle \boldsymbol{v}|\boldsymbol{u}\rangle\).

Note that, given a cell pair \((i,j)\) or a spike train pair \((\boldsymbol{u},\boldsymbol{v})\), \(c_{ij}\) does not depend on spike times and \(\langle \boldsymbol{u}|\boldsymbol{v}\rangle\) does not depend on cell labeling.

With this notation, we can define the multi-unit spike train distance as

\[\label{eq:distance_as_intprod} \left\Vert |\mathcal{U}\rangle - |\mathcal{V}\rangle \right\Vert^2 = \langle \mathcal{U}|\mathcal{U}\rangle + \langle \mathcal{V}|\mathcal{V}\rangle - 2 \langle \mathcal{U}|\mathcal{V}\rangle\]

where the multi-unit spike train inner product \(\langle \mathcal{V}|\mathcal{U}\rangle\) between \(\mathcal{U}\) and \(\mathcal{V}\) is just the natural bilinear operation induced on \(\mathcal{S}\) by the tensor product structure:

\[\begin{split}\label{eq:multiunit_intprod} \begin{split} \langle \mathcal{V}|\mathcal{U}\rangle &= \sum_{i,j=1}^C \langle i|j \rangle\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle = \sum_{i,j=1}^C c_{ij}\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle\\ &= \sum_{i=1}^C\left[ c_{ii} \langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c_{ij} \left(\sum_{j<i}\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle + \sum_{j>i}\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle \right) \right] \end{split}\end{split}\]

But \(c_{ij}=c_{ji}\) and \(\langle \boldsymbol{v}|\boldsymbol{u}\rangle=\langle \boldsymbol{u}|\boldsymbol{v}\rangle\), so

\[\begin{split}\begin{split} \sum_{i=1}^C\sum_{j<i}c_{ij}\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle &= \sum_{j=1}^C\sum_{i<j}c_{ji}\langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle = \sum_{i=1}^C\sum_{j>i}c_{ji}\langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle = \sum_{i=1}^C\sum_{j>i}c_{ij}\langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle\\ &=\sum_{i=1}^C\sum_{j>i}c_{ij}\langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle \end{split}\end{split}\]

and

\[\begin{split}\label{eq:multiprod_j_ge_i} \langle \mathcal{V}|\mathcal{U}\rangle = \sum_{i=1}^C\left[ c_{ii} \langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c_{ij} \sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle \right) \right]\end{split}\]

Now, normally we are interested in the particular case where \(c_{ij}\) is the same for all pair of distinct cells:

\[\begin{split}c_{ij} = \begin{cases} 1 & \textrm{if } i=j\\ c & \textrm{if } i\neq j \end{cases}\end{split}\]

and under this assumption we can write

\[\begin{split}\label{eq:multiprod_constant_c} \langle \mathcal{V}|\mathcal{U}\rangle = \sum_{i=1}^C\left[\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle \right) \right]\end{split}\]

and

\[\begin{split}\begin{split} \left\Vert |\mathcal{U}\rangle - |\mathcal{V}\rangle \right\Vert^2 = \sum_{i=1}^C\Bigg\{&\langle \boldsymbol{u}^i|\boldsymbol{u}^i \rangle + c\sum_{j>i}\left(\langle \boldsymbol{u}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{u}^j \rangle \right) +\\ +&\langle \boldsymbol{v}^i|\boldsymbol{v}^i \rangle + c\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{v}^j \rangle + \langle \boldsymbol{v}^i|\boldsymbol{v}^j \rangle \right) + \\ -2&\left[\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle \right)\right] \Bigg\} \end{split}\end{split}\]

Rearranging the terms

\[\begin{split}\begin{gathered} \label{eq:multidist_constant_c} \left\Vert |\mathcal{U}\rangle - |\mathcal{V}\rangle \right\Vert^2 = \sum_{i=1}^C\Bigg[\langle \boldsymbol{u}^i|\boldsymbol{u}^i \rangle + \langle \boldsymbol{v}^i|\boldsymbol{v}^i \rangle - 2 \langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + \\ + 2c\sum_{j>i}\left(\langle \boldsymbol{u}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{v}^i|\boldsymbol{v}^j \rangle -\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle - \langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle\right)\Bigg]\end{gathered}\end{split}\]

Van Rossum-like metrics

In Van Rossum-like metrics, the feature function and the single-unit kernel are, for \(\tau\neq 0\),

\[\begin{split}\begin{gathered} \phi^{\textrm{VR}}_{\tau}(t) = \sqrt{\frac{2}{\tau}}\cdot e^{-t/\tau}\theta(t) \\ \mathcal{K}^{\textrm{VR}}_{\tau}(t_1,t_2) = \begin{cases} 1 & \textrm{if } t_1=t_2\\ e^{-\left\vert t_1-t_2 \right\vert/\tau} & \textrm{if } t_1\neq t_2 \end{cases}\end{gathered}\end{split}\]

where \(\theta\) is the Heaviside step function (with \(\theta(0)=1\)), and we have chosen to normalise \(\phi^{\textrm{VR}}_{\tau}\) so that

\[\left\Vert \phi^{\textrm{VR}}_{\tau} \right\Vert_2 = \sqrt{\int\textrm{d }\!t \left[\phi^{\textrm{VR}}_{\tau}(t)\right]^2} = 1 \quad.\]

In the \(\tau\rightarrow 0\) limit,

\[\begin{split}\begin{gathered} \phi^{\textrm{VR}}_{0}(t) = \delta(t)\\ \mathcal{K}^{\textrm{VR}}_{0}(t_1,t_2) = \begin{cases} 1 & \textrm{if } t_1=t_2\\ 0 & \textrm{if } t_1\neq t_2 \end{cases}\end{gathered}\end{split}\]

In particular, the single-unit inner product now becomes

\[\langle \boldsymbol{u}|\boldsymbol{v}\rangle = \sum_{n=1}^N\sum_{m=1}^M \mathcal{K}^{\textrm{VR}}(u_n,v_m) = \sum_{n=1}^N\sum_{m=1}^M e^{-\left\vert u_n-v_m \right\vert/\tau}\]

Markage formulas

For a spike train \(\boldsymbol{u}\) of length \(N\) and a time \(t\) we define the index \(\tilde{N}\left( \boldsymbol{u}, t\right)\)

\[\begin{split}\tilde{N}\left( \boldsymbol{u}, t\right) \doteq \max\{n | u_n < t\}\end{split}\]

which we can use to re-write \(\langle \boldsymbol{u}|\boldsymbol{u}\rangle\) without the absolute values:

\[\begin{split}\begin{split} \langle \boldsymbol{u}|\boldsymbol{u}\rangle&= \sum_{n=1}^N \left( \sum_{m|v_m<u_n}e^{-(u_n-v_m)/\tau} + \sum_{m|v_m>u_n}e^{-(v_m-u_n)/\tau} + \sum_{m=1}^M\delta\left(u_n,v_m\right) \right)\\ &= \sum_{n=1}^N \left( \sum_{m|v_m<u_n}e^{-(u_n-v_m)/\tau} + \sum_{m|u_m<v_n}e^{-(v_n-u_m)/\tau} + \sum_{m=1}^M\delta\left(u_n,v_m\right) \right)\\ &= \sum_{n=1}^N \left[ \sum_{m=1}^{\tilde{N}\left( \boldsymbol{v}, u_n\right)}e^{-(u_n-v_m)/\tau} + \sum_{m=1}^{\tilde{N}\left( \boldsymbol{u}, v_n\right)}e^{-(v_n-u_m)/\tau} + \delta\left(u_n,v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)+1}\right) \right]\\ &= \sum_{n=1}^N \Bigg[ e^{-(u_n-v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)})/\tau} \sum_{m=1}^{\tilde{N}\left( \boldsymbol{v}, u_n\right)}e^{-(v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)}-v_m)/\tau} + \\ &\phantom{ = \sum_{n=1}^N } + e^{-(v_n-u_{\tilde{N}\left( \boldsymbol{u}, v_n\right)})/\tau} \sum_{m=1}^{\tilde{N}\left( \boldsymbol{u}, v_n\right)}e^{-(u_{\tilde{N}\left( \boldsymbol{u}, v_n\right)}-u_m)/\tau} + \\ &\phantom{ = \sum_{n=1}^N } + \delta\left(u_n,v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)+1}\right) \Bigg]\\ \end{split}\end{split}\]

For a spike train \(\boldsymbol{u}\) of length \(N\), we also define the the markage vector \(\boldsymbol{m}\), with the same length as \(\boldsymbol{u}\), through the following recursive assignment:

\[\begin{split}\begin{aligned} m_1(\boldsymbol{u}) &\doteq 0 \\ m_n(\boldsymbol{u}) &\doteq \left(m_{n-1} + 1\right) e^{-(u_n - u_{n-1})/\tau} \quad \forall n \in \{2,\ldots,N\}\label{eq:markage_definition}\end{aligned}\end{split}\]

It is easy to see that

\[\begin{split}\begin{split} m_n(\boldsymbol{u}) &= \sum_{k=1}^{n-1}e^{-(u_n - u_k)/\tau} = \left(\sum_{k=1}^{n}e^{-(u_n - u_k)/\tau}\right) - e^{-(u_n - u_n)/\tau}\\ &= \sum_{k=1}^{n}e^{-(u_n - u_k)/\tau} - 1 \end{split}\end{split}\]

and in particular

\[\label{eq:markage_sum} \sum_{n=1}^{\tilde{N}\left( \boldsymbol{u}, t\right)}e^{-(u_{\tilde{N}\left( \boldsymbol{u}, t\right)}-u_n)/\tau} = 1 + m_{\tilde{N}\left( \boldsymbol{u}, t\right)}(\boldsymbol{u})\]

With this definition, we get

\[\begin{split}\label{eq:singleunit_intprod_markage} \begin{split} \langle \boldsymbol{u}|\boldsymbol{v}\rangle &= \sum_{n=1}^N \Bigg[ e^{-(u_n-v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)})/\tau} \left(1 + m_{\tilde{N}\left( \boldsymbol{v}, u_n\right)}(\boldsymbol{v})\right) + \\ &\phantom{= \sum_{n=1}^N} + e^{-(v_n-u_{\tilde{N}\left( \boldsymbol{u}, v_n\right)})/\tau} \left(1 + m_{\tilde{N}\left( \boldsymbol{u}, v_n\right)}(\boldsymbol{u})\right) + \\ &\phantom{= \sum_{n=1}^N} + \delta\left(u_n,v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)+1}\right) \Bigg] \end{split}\end{split}\]

Finally, note that because of the definition of the markage vector

\[e^{-(u_n-u_{\tilde{N}\left( \boldsymbol{u}, u_n\right)})/\tau} \left(1 + m_{\tilde{N}\left( \boldsymbol{u}, u_n\right)}(\boldsymbol{u})\right) = e^{-(u_n-u_{n-1})/\tau}\left(1+m(\boldsymbol{u})\right) = m_{n}(\boldsymbol{u})\]

so that in particular

\[\begin{split}\label{eq:singleunit_squarenorm_markage} \begin{split} \langle \boldsymbol{u}|\boldsymbol{u}\rangle &= \sum_{n=1}^N \left(1 + 2m_{n}(\boldsymbol{u})\right) \end{split}\end{split}\]

A formula for the efficient computation of multiunit Van Rossum spike train metrics, used by pymuvr, can then be obtained by opportunely substituting these expressions for the single-unit scalar products in the definition of the multiunit distance.