Large amplitude spikes classified incorrectly into multiple clusters

[samamckenzie]

This issue is best illustrated with a picture (see attached)

Kilosort originally gave me two units (yellow and green in the attached plot) but I noticed that there was unphysiologically high 0ms correlation. I then noticed that the green unit showed two distinct clusters which, when manually split revealed the source of the high correlation (lavender cluster). Notice in the attached cross correlations that the lavender and yellow show many 0ms correlations and share a refractory period: they are the same spikes. Lavender and green do not share a refractory period: they must originate from different neurons.

This issue happens every time there is a large amplitude unit (or any large amplitude noise, actually) and is easily spotted by identifying: unusually high 0ms correlations without any jitter, bimodal waveform clusters that, when split, reveal an absence of a common refractory period. This manual stage is, in my opinion, essential for the software to be unusable for our purposes, and easily adds double the manual sorting time as compared to the klustakwik pipeline. Failing to inspect and correct for this issue, unfortunately, takes the best signal (large amplitude spikes) and distributes it across many units, making it the worst signal for looking at the correlations between single units.

I understand the root of the problem: variance and mean scale together, so large amplitude spikes are the most variable meaning that the residuals will be highest when looking for template (mis)matches. Perhaps it is possible to adjust the inclusion threshold according to spike amplitude?

thank you sam double_spike.pdf

[mnpompili]

Hi Sam,

we are encountering the same problem here at the Zugaro lab. Many large amplitude spikes assigned to multiple clusters (especially for hippocampal recordings of pyramidal cells). Did you end up finding an explanation/solution for this?

thanks, Marco

[samamckenzie]

Dear Marco

Explanation: Waveform changes over burst. Mean and variance scale together. Mean drops with square of distance. So nearby cells will have a very poor template that will match to waveforms that are smaller than the first in the burst. When substracted from the early spikes there will be energy left over that is then fit to a smaller unit with a waveform on the same channel(s). The issue is always between the best signal and the worst

Solution: You need to look at every pairwise CCG for every unit on the same shank. When you see the 0lag no jitter synchrony, you can bet this is the problem. Two solutions. 1) recluster these pairs with klustakwik. 2) manually find the odd large spikes on the small unit. Both require you to recalculate the PCs in a common space. I cannot stand phy. since 1) it isn't naturally organized by shank and 2) the PCs are local and therefore plotting two units together is meaningless. Of course it is open source and both can be fixed, but I am still using klusters to clean up the units. we have conversion scripts if you need them

sam

On Wed, Nov 21, 2018 at 5:38 AM mnpompili notifications@github.com wrote:

Hi Sam,

we are encountering the same problem here at the Zugaro lab. Many large amplitude spikes assigned to multiple clusters (especially for hippocampal recordings of pyramidal cells). Did you end up finding an explanation/solution for this?

thanks, Marco

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/cortex-lab/KiloSort/issues/47#issuecomment-440617100, or mute the thread https://github.com/notifications/unsubscribe-auth/ANnfgMJDmKgzVDjPLF1xrP8rrGO96vZiks5uxS0xgaJpZM4LUOlt .

[mnpompili]

Hi Sam,

we also wrote conversation scripts to keep using klusters and neuroscope. It looks like either us or you could have avoided it... We have this problem of simultaneous spikes a lot (see below, e.g. clusters 4 and 6 have a peak at 0ms with all other units) It looks that the problem is reduced when forcing kilosort to split more its clusters; I'm running some tests now but I also have the problem of not enough memory on the GPU if I asks for too many clusters.

Another option would be to not allow kilosort to detect simultaneous spikes. We could not find an option to disable it but maybe @marius10p would know whether this is possible?

[samamckenzie]

Yes to more clusters. You can also force the classification closer to the mean. You don't want to eliminate the overlap, that will defeat the purpose of kilosort.

Re GPU memory, you have this problem even when you decrease the GPU memory in the param file? We have never had this problem even with 256 channels. What model GPU. We can send a param file if you need one.

All the best Sam

On Nov 21, 2018 8:36 AM, "mnpompili" notifications@github.com wrote:

Hi Sam,

we also wrote conversation scripts to keep using klusters and neuroscope. It looks like either us or you could have avoided it... We have this problem of simultaneous spikes a lot (see below, e.g. clusters 4 and 6 have a peak at 0ms with all other units) [image: image] https://user-images.githubusercontent.com/45228967/48844354-e4405280-ed99-11e8-8142-d27b95a79614.png It looks that the problem is reduced when forcing kilosort to split more its clusters; I'm running some tests now but I also have the problem of not enough memory on the GPU if I asks for too many clusters.

Another option would be to not allow kilosort to detect simultaneous spikes. We could not find an option to disable it but maybe @marius10p https://github.com/marius10p would know whether this is possible?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/cortex-lab/KiloSort/issues/47#issuecomment-440664071, or mute the thread https://github.com/notifications/unsubscribe-auth/ANnfgC4I7jArfiypnUhhCX6MQJoEWm26ks5uxVbHgaJpZM4LUOlt .

[mnpompili]

I know that one of the main features of kilosort is the detection of overlapping spikes but we are willing to drop them to solve this issue of big spikes contaminating all clusters. The main reason why I am using kilosort is because it's faster. KlustaKwik on my data (1 session = 9h of recording, 256ch) takes an inifinite time (I dropped it after more than 1 month of clustering for one session).

About the GPU we use the 1080ti. By GPU memory parameter do you mean the ops.NT? We use 64*1024+ ops.ntbuff.

I am happy to take your parameter file to compare it with ours (below)

All the best, Marco

`function ops = KilosortOptions(datFile,channelMapFile) % Load options for Kilosort

if any(strfind(datFile,'/')), slash = '/'; else slash = '\'; end [directory,~] = fileparts(datFile); directory = [directory slash];

useGPU = 1; % index of GPU to use. If == 0, uses CPU

ops.fpath = directory; ops.GPU = useGPU; % whether to run this code on an Nvidia GPU (much faster, mexGPUall first): index of GPU to use ops.parfor = 1; % whether to use parfor to accelerate some parts of the algorithm ops.verbose = 1; % whether to print command line progress ops.showfigures = 0; % whether to plot figures during optimization

ops.datatype = 'dat'; % binary ('dat', 'bin') or 'openEphys' ops.fbinary = datFile; % will be created for 'openEphys' ops.fproc = fullfile(directory, 'temp_whitened_signal.dat'); % residual from RAM of preprocessed data ops.root = directory; % 'openEphys' only: where raw files are % define the channel map as a filename (string) or simply an array ops.chanMap = channelMapFile; % make this file using createChannelMapFile.m % ops.chanMap = 1:ops.Nchan; % treated as linear probe if unavailable chanMap file

ops.Nfilt = 960; % number of clusters to use (2-4 times more than Nchan, should be a multiple of 32) (GM: MAX WORKING: 960 // 1152 crashes GPU) ops.nNeighPC = []; % visualization only (Phy): number of channnels to mask the PCs, leave empty to skip (12) ops.nNeigh = 32; % visualization only (Phy): number of neighboring templates to retain projections of (16)

% options for channel whitening ops.whitening = 'full'; % type of whitening (default 'full', for 'noSpikes' set options for spike detection below) ops.nSkipCov = 1; % compute whitening matrix from every N-th batch (1) ops.whiteningRange = 32; % how many channels to whiten together (Inf for whole probe whitening, should be fine if Nchan<=32) // should try 8

ops.criterionNoiseChannels = 0.1; % fraction of "noise" templates allowed to span all channel groups (see createChannelMapFile for more info).

% other options for controlling the model and optimization ops.Nrank = 3; % matrix rank of spike template model (3) ops.nfullpasses = 6; % number of complete passes through data during optimization (6) ops.maxFR = 20000; % maximum number of spikes to extract per batch (20000) ops.fshigh = 300; % frequency for high pass filtering % ops.fslow = 2000; % frequency for low pass filtering (optional) ops.ntbuff = 64; % samples of symmetrical buffer for whitening and spike detection ops.scaleproc = 200; % int16 scaling of whitened data ops.NT = 641024+ ops.ntbuff;% this is the batch size (try decreasing if out of memory) 641024+ ops.ntbuff = 500MB of VRAM % for GPU should be multiple of 32 + ntbuff

% the following options can improve/deteriorate results. % when multiple values are provided for an option, the first two are beginning and ending anneal values, % the third is the value used in the final pass. ops.Th = [3 6 6]; % threshold for detecting spikes on template-filtered data ([6 12 12]) ([5 10 10]) ([4 8 8]) // Penalty for spike detection. The higher it is; the less spikes detected ops.lam = [10 30 30]; % large means amplitudes are forced around the mean ([10 30 30]) // The higher these values are, the more Kilosort ensures spikes within a cluster have close amplitudes ops.nt0 = 5; % Number of samples to take around a spike ops.nannealpasses = 4; % should be less than nfullpasses (4) ops.momentum = 1./[20 800]; % start with high momentum and anneal (1./[20 1000]) ops.shuffle_clusters = 1; % allow merges and splits during optimization (1) ops.mergeT = .1; % upper threshold for merging (.1) ops.splitT = .1; % lower threshold for splitting (.1)

% options for initializing spikes from data ops.initialize = 'no'; %'fromData' or 'no' ops.spkTh = 2; % spike threshold in standard deviations (4) ---------------------------------- ops.loc_range = [3 1]; % ranges to detect peaks; plus/minus in time and channel ([3 1]) ops.long_range = [30 6]; % ranges to detect isolated peaks ([30 6]) ops.maskMaxChannels = 5; % how many channels to mask up/down ([5]) ops.crit = .65; % upper criterion for discarding spike repeates (0.65) ops.nFiltMax = 10000; % maximum "unique" spikes to consider (10000)

% load predefined principal components (visualization only (Phy): used for features) dd = load('PCspikes2.mat'); % you might want to recompute this from your own data ops.wPCA = dd.Wi(:,1:7); % PCs

% options for posthoc merges (under construction) ops.fracse = 0.1; % binning step along discriminant axis for posthoc merges (in units of sd) ops.epu = Inf;

ops.ForceMaxRAMforDat = 32e9; % maximum RAM the algorithm will try to use; on Windows it will autodetect.`

[brendonw1]

I had an idea in relation to the big spikes problem. Basically this is due to a desire to catch burst spikes, right? But of course burst spikes are smaller than the modal spike... the big spikes are of course bigger than the mode. So one idea is to find the mode, and allow in putative burst spikes smaller than that mode but then use some cutoff to exclude large spikes above that mode/peak of distribution. Of course you have to allow some spikes bigger than the peak, but there may be some disciplined way of looking at the right side of the distribution and finding a cutoff that for the most part allows in real spikes and excludes huge ones. This also fits with theory.

By the way I also think we may want to actually go back after all is done (even klusters/phy) and have an algo specifically hunt down burst spikes by specifically looking for reduced-amplitude versions of initial spikes at a few ms delay.

Finally, I modified klusters to allow it to work with >24hr data sets (wasn't displaying time axis right so couldn't look at stability)... if anyone wants that version let me know.

Brendon

[petersenpeter]

In my experience, I believe the biggest issue is with artifacts. With artifacts, you can have other units "templated" on top of the noise. Between similar sized units I rarely see it. mnpompili, I believe it could be related to your very low threshold for spike detection (ops.Th = [3 6 6]). From my experience having it below [6 10 10] is not necessary/good as you end up with a ton of clusters that cannot be cleaned anyhow. These events where Kilosort adds the same spike time stamps to multiple clusters happens at large artifacts, and you should toss the noise artifacts as these are the spurious timestamps contaminating your CCGs. Further, if you experience the problem with similar sized waveforms, I believe it has to do with how you are doing your cluster cutting. Lets say that Kilosort creates two templates that have matching timestamps, meaning that two units' spiking overlaps in time, if you separate the first unit into two separate units, one of these separated units will have matching timestamps with the secondary original cluster. So when cluster cutting, you should often only keep the main cluster, and discard the other spikes as they are already captured with another template.

[mnpompili]

Hi Peter We had to lower the threshold because with [5 10 10] and also [4 8 8] a lot of spikes were not extracted. With [3 6 6] we actually do not get too many clusters nor artifact spikes: on a given electrodes bundle (4 to 8 channels) on average we end up with 20 good clusters cleaned out of 40 provided by kilosort. I have to specify that 1) we have a script that get rid of artifact clusters before we actually open the clu file in klusters so it makes cleaning way easier 2) before running kilosort we detect large artifact events like the animal bumping its head etc.. and we set the signal during these events to zero to avoid these very large event contaminating our clusters.

My problem is rather that I would like to split even more the clusters I obtain because they are contaminated by spikes belonging to other cells. However I cannot do it because my GPU runs out of memory if I ask for more clusters (see above).

Therefore we do not think that these zero lag Xcorr come from artifact spikes. And I do not think it's a problem coming from my treshold. Below you see the same clusters extracted with a treshold at [5 10 10] and the zero lag correlations are still there

[brendonw1]

What is the script for cleaning noise spikes ahead of time? Something you'd want to share? I think kilosort is great but is a bit lightly supported considering Marius et al are genuinely busy. So it's be great to come together as a community as much as possible... git repos etc. (I have some ideas myself but getting recordings going in the lab before focusing on sorting details)

[mnpompili]

we are of course planning to share our functions including these cleaning functions but we are still testing them and given the results above we are still not completely happy with these. Also the code is still a complete mess. But we will definitely share this once we finish our tests and tidy the code up.

[brendonw1]

Totally understood of course... thanks!

[mnpompili]

To come back to the original issue, our current plan is to write a script to get rid of these. This because we are interested in cells that fires together and this clustering behavior will create a lot of perfect co-firing between clusters that would be mostly artifactual. Our script would look for each spike associated to multiple clusters and assign it to the cluster with the most similar average waveform. It seems fairer to us to decide that a given peak belongs only to one cell.

Is anybody aware of some setting in the Kilosort native scripts to avoid it detecting simultaneous spikes? Such a setting would spare us the time to develop the script described above....

[brendonw1]

The only think I know of that would minimiize almost-simultaneous spikes is to reduce the temporal duration of the template... but that would only handle double spikes more distant in time.

We'd also thought of a script to reassign the double spikes to at least one of the templates... as I think you're implying you could, if sufficiently confident also assign a lagged spike to a second cluster. It'd be great for the community to have that function.

Maybe we could make a community repository of add-ons to Kilosort?... probably mostly matlab based since kilosort is matlab-based?

On Tue, Nov 27, 2018 at 5:08 AM Marco Pompili notifications@github.com wrote:

To come back to the original issue, our current plan is to write a script to get rid of these. This because we are interested in cells that fires together and this clustering behavior will create a lot of perfect co-firing between clusters that would be mostly artifactual. Our script would look for each spike associated to multiple clusters and assign it to the cluster with the most similar average waveform. It seems fairer to us to decide that a given peak belongs only to one cell.

Is anybody aware of some setting in the Kilosort native scripts to avoid it detecting simultaneous spikes? Such a setting would spare us the time to develop the script described above....

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/cortex-lab/KiloSort/issues/47#issuecomment-442001374, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTUF491oHERobK38AvN48Sp5x-xNOks5uzQ8bgaJpZM4LUOlt .