# 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.