By: darioringach
M: Analysis, Uncategorized

Spikefinder: the Vanilla algorithm

Over the last few months, the Spikefinder challenge has provided a playing ground for colleagues to offer ideas about how to estimate the spiking of individual neurons based on the measured activity of fluorescent calcium indicators.

The challenge was for people to come up with strategies that beat the performance of state-of-the-art algorithms, STM & oopsi. A good number of algorithms were able to achieve this in short time, including one I submitted, termed Vanilla.

The best performing algorithms seem to have relied on modern machine learning methods. Vanilla is nothing more than a linear filter followed by a static-nonlinearity y(t) = \phi ( h(t) * x(t) )  — thus the name.

The filter h(t) is a linear combination of an even filter, estimating the mean of the signal at time t, and an odd filter, estimating the derivative of the signal at time t.

The even filter is a Gaussian, h_{even} = A \exp ( -t^2 / 2 \sigma^2 ), and the odd filter is the derivative of a Gaussian h_{odd} = B t \exp (-t^2 / 2 \sigma^2).  The constants A and B are such that the norm of the filters is normalized to one, \| A \| = \| B \| = 1.  These two filters are linearly combined while keeping the norm of resulting filter equal to one, h(t) = \cos \alpha \: h_{even}(t) + \sin \alpha \: h_{odd}(t).

The output nonlinearity is a rectifier to a power, \phi ( x ) = (x- \theta)^\beta if x>\theta, and zero otherwise.

The model has only 4 parameters, \{\sigma, \alpha, \theta, \beta \}. The amount of smoothing of the signal is controlled by \sigma, the shape of the filter is controlled by \alpha, and the threshold \theta and power \beta determine the shape of the nonlinearity.

The model is fit by finding the optimal values of \{\sigma, \alpha, \theta, \beta \} that maximize the correlation between its output y(t) and the recorded spiking of the neuron.  I used Matlab’s fminsearch() to perform this optimization, which was typically finished in about 60 sec or less for most datasets.

The only pre-processing done was a z-scoring of the raw signals.  In one dataset (dataset #5, GCaMP6s in V1), we allowed for an extra-delay parameter between the signal and the prediction.

I was surprised this entry did relatively well, as the algorithm it is basically a version of STM. I think the particular shape of the output nonlinearity (rectifier+power vs exponential), the constrain imposed on the shape of the filters, and the resulting small number of parameters, paid a role in Vanilla doing better overall.

The top algorithms reached an absolute performance of about 0.47 and it seems unlikely this performance can be improved by a lot. This seems to highlight the limitations of the current families of calcium indicators in yielding precise spiking information — so there is plenty of opportunity to improve them.

It is interesting that despite its simplicity, the relative performance of Vanilla, with a correlation coefficient of 0.428 was not dramatically inferior to that of the top performing, deep network models, with all its bells and whistles, which landed at 0.464.  So, one must pay due respects to deep networks, but I was honestly expecting Vanilla to be completely blown out of the water in terms of performance, and I don’t think it was.

Finally, Vanilla is a non-casual algorithm, as both past and future samples are used to predict the response at time t. In some situations, however, when we are trying close the loop as fast as possible by controlling a stimulus based on neural activity itself, we need algorithms that are casual and can provide a fast, online estimate of spiking activity. I wonder if any of the submissions are causal algorithms and, if not, what would be the best performance the methods can attain if we allow them to provide estimates of spiking based only on past samples.