Fixing upsampled spike count and alpha

[kaghi]

What exactly was happening?
During the upsampling step spikes can distributed into k-bins with some centroid. FISTA performs a grid search using a Box[0,1] and distributes a single spike's upsampled response across adjacent bins to the centroid rather than keep the upsampled spikes constrained to a single bin. Once this step occurs we then binarize every bin and threshold as its own spike. Downsample_binary then sums them into one frame at base sampling frequency which turns one spike into k spikes.

Why does this matter?
K·s is the reconvolved calcium train. Even after upsampling by 10x the shifted kernel is basically identical to the unshifted kernel. So for the reconvolved trace this looks like:

K·s = (s[i] + s[i+1] + s[i+2] + …) · (whatever your kernel shape is). This K·s is the mass we will integrate.

The issue is that the convolution can't really see how this mass is split across the bins, it can only obtain the sum. In other words, you could technically have two spike trains with the same mass but entirely different distributions (eg. a flat set of bins vs one tall bin) and you'd get the same K·s. This results in an identical reconstruction error (PVE) that is flat regardless of the threshold -> count outcome you choose IF and only if the mass is preserved. However, the issue is that despite preserved mass you get wildly different s_counts. The timing is correct and is not an issue with proper inference. This paper by Duval-Peyre basically predicts this particular phenomenon and recommends working from the aggregate mass: https://hal.science/hal-01135200v2/document.

How do we go about fixing this?

Instead of counting bins above a threshold, we can integrate under the mass curve (the spread out upsampled continuous spike mass across bins) and then we can convert that to a spike count. For each run we send a single synthetic spike and measure the mass it generates when pushed through the same deconvolution. This is our calibration mass (calibration_mass). The s_count is then very simply the run_mass/calibration_mass and rounded. We send a single spike every time we push data through the FISTA solver so both the synthetic spike and real data are processed exactly the same.

To make sure sure we are integrating proper spike mass and not noise we have to establish some fixed minimum floor before integration (0.5 / upsample_factor).min(0.15).

Once the count is right we can then refit the alpha. It looks like it works pretty well in recovering alpha.