Adapt template metrics to work with positive spiking neurons (peak_sign='pos')

[ecobost]

Template metrics are currently designed to only work for negative spiking neurons (large negative trough followed by small positive peak) but many units spike higher in the positive side (very common in NHP). I adapted here those metrics to deal properly with positive spiking units when peak_sign='pos' (or 'both'). When peak_sign='neg' (default), results are the same.

For reference, here are the templates for a positive spiking unit (each line is a channel).

Currently, the cyan line will be picked as the template, trough will be at ~40 (the dip) and peak will be the largest value after trough (so probably that small peak at 60 or the one at 75). That will produce non-sensical metrics (in this case, peak-to-through distance will be overestimated, halfwidth will be overestimated, peak-to-trough valley might be understimated). This is true for all positive spiking units. I propose that when the unit spikes upwards (as declared by peak_sign='pos'), the positive peak should be the "trough" and the negative dip post-positive spike should be the "peak", which will result in metrics that make sense for this units and is consistent with how peak_sign='pos' is dealt with in other parts of spikeinterface. I don't love the overriding of the name "trough" and "peak" (maybe there are more general terms for that but I didn't wanna change much of the code and the names make sense for the default, more common, peak_sign='neg')

Most metrics work fine whether the spike is positive or negative. A few needed some tweaks; I describe what I changed here (per metric):

• halfwidth: halfwidth: currently implemented as find the first index before through and last index after through where template is less than threshold; in this PR it is implemented as find the last index before trough and first index after trough where the threshold is crossed (regardless of direction). This computes the halfwidth of the spike ("trough") whether it points downwards or upwards. The slight logic change in how indices are picked is intentional, it deals better with the edge case where another spike in the same template also crosses the threshold (see image)*
  
• repolarization_slope: current implementation finds the first index after trough that is greater than zero; this PR changes it to find the first index after trough where the zero is crossed (which works for negative and positive spikes). if using peak_sign='both', the slope will then be positive for negative/downward spiking units and negative for upwards spiking units, which is nice.
• velocity_fits: I invert the templates if using peak_sign='pos' so it fits the speed of the positive peaks rather than the negative peaks.

Note: in halfwidth, there is another small change that will affect peak_sign='neg' data. Currently, if a template was [0, -1, 0, 0, 0], code will say the halfwidth was 2, if a template was [0, -1, -1, 0, 0] it will say the half_width is 3 and so on. So there was a one-off error. I fixed that so now if x indices are above the threshold, the halfwidth is x (rather than x + 1), but then most new data will differ by one from previous data (e.g., 30 microseconds shorter halfwidth for neuropixels)

[alejoe91]

HI @ecobost

We are entirely refactoring the find peaks and template metrics in this other PR #4306

I think that the changes you propose here are included (and expanded) there. Can you take a look and give it a try?

[ecobost]

Hi @alejoe91, #4306 adds extra metrics but it still assumes all units spike downwards: it still assigns the biggest negative peak ("through") as the main peak and the positive peaks as the secondary ones (out of which, halfwidth, peak_to_through_duration, repolarization_slope, recovery_slope, waveform_widths, waveform_ratios and velocity_fits are computed). This gives results that don't make sense for positive-spiking units (e.g., halfwidth of the largest negative peak rather than the positive peak, repolarization slope computed from the largest negative peak to zero rather than from the largest positve peak to zero, waveform ratios between the trough and positive peaks to either side rather than the positive peak and troughs to either side and so on).

My idea was to rephrase template metrics as computed from a main peak (could be positive or negative) and secondary/adjacent peaks (which are the opposite direction of the main peak). For upwards spiking units, this is more informative. And peak_sign would control the direction of the main peak. This is what this PR is doing. Something similar can be done in #4306.

This will be very useful, more for NHP where positive spiking units are common and templates look like the flipped version of negative spiking units. And once bombcell gets incorporated, it will make classification work (hopefully) for this positve-spiking units too (which I presume currently will all be classified poorly as the metrics are unreliable).

[alejoe91]

Hi @alejoe91, #4306 adds extra metrics but it still assumes all units spike downwards: it still assigns the biggest negative peak ("through") as the main peak and the positive peaks as the secondary ones (out of which, halfwidth, peak_to_through_duration, repolarization_slope, recovery_slope, waveform_widths, waveform_ratios and velocity_fits are computed). This gives results that don't make sense for positive-spiking units (e.g., halfwidth of the largest negative peak rather than the positive peak, repolarization slope computed from the largest negative peak to zero rather than from the largest positve peak to zero, waveform ratios between the trough and positive peaks to either side rather than the positive peak and troughs to either side and so on).

My idea was to rephrase template metrics as computed from a main peak (could be positive or negative) and secondary/adjacent peaks (which are the opposite direction of the main peak). For upwards spiking units, this is more informative. And peak_sign would control the direction of the main peak. This is what this PR is doing. Something similar can be done in #4306.

This will be very useful, more for NHP where positive spiking units are common and templates look like the flipped version of negative spiking units. And once bombcell gets incorporated, it will make classification work (hopefully) for this positve-spiking units too (which I presume currently will all be classified poorly as the metrics are unreliable).

@ecobost I think you're right! @Julie-Fabre in bombcell and in #4306 you always compute the halfwidth and peak to trough duration from the main trough. Is this intented? Or should we comput the halfwidth on the main peak (which could be positive) and also use either pre or post positive peaks for peak to trough duration?