Unit Quality Metrics

Overview

This Jupyter notebook will provide a detailed explanation of the unit quality metrics included the various Neuropixels datasets. These quality metrics enable you to identify units that are well separated, single unit data. It’s important to pay attention to quality metrics, because failing to apply them correctly could lead to invalid scientific conclusions, or could end up hiding potentially useful data.

To help you avoid these pitfalls, this tutorial will explore how these metrics are calculated, how they can be biased, and how they should be applied to specific use cases. It’s important to keep in mind that none of these metrics are perfect, and that the use of unit quality metrics for filtering ephys data is still an evolving area of research. More work is required in order to establish general-purpose best practices and standards in this domain.

Why do we need quality metrics?

For a long time, converting continuous voltage traces to sorted spike times was one of the “dark arts” of neuroscience. Spikes were typically sorted by hand-drawing boundaries around clouds of dots, using heuristics learned from other lab members. The quality of the resulting clusters could be validated by looking at metrics such as ISI violations or isolation distance, but there were no standards governing how these metrics informed which units to include for further analysis.

Recent in advances in neural recording devices, such as Neuropixels, have made it practically impossible to sort spikes by hand. Fortunately, we now have access to powerful algorithms that use GPUs to sort spikes in approximately the same amount of time it took to record the data. All of the Allen Institute Neuropixels data has been sorted with Kilosort, a template-matching algorithm developed by Marius Pachitariu at HHMI Janelia Research Campus.

For Neuropixels recordings with minimal electrode drift, Kilosort2 performs well enough that further manual curation is not necessary. Unlike the original version of Kilosort, which required a manual merging step, Kilosort2 attempts to merge units automatically. Sometimes it over-merges, leading to units that clearly combine spikes from multiple cells. But in the majority of cases, Kilosort2 makes merging decisions as well as a human would, and does so in a way that is highly reproducible.

Because there is no “ground truth” information available in these datasets, any sorting algorithm is bound to make mistakes. Quality metrics allow us to understand the types of mistakes that are occurring, and obtain an estimate of their severity. Some common errors that can be identified by quality metrics include:

• Assigning spikes from multiple neurons to the same cluster

• Missing spikes from neurons with waveform amplitude near the spike detection threshold

• Failing to track neurons with waveforms that change as a result of electrode drift

These mistakes can occur even in units that appear to be extremely well isolated. It’s misleading to conceive of units as existing in two distinct categories, one with perfectly clean “single units” and one with impure “multiunits.” Instead, there’s a gradient of qualities, with mostly complete, uncontaminated units at one end, and incomplete, highly contaminated units at the other.

Despite the fact that there’s not a clear division between single-unit and multi-unit activity, we still have to make a binary decision in every analysis we carry out: should this unit be included or not? Ideally this decision should be based on objective metrics that will not bias the end results. Our usual default qc uses three quality metrics, isi_violations, amplitude_cutoff, and presence_ratio, to filter out units that are likely to be highly contaminated or missing lots of spikes. However, the default values of these filters may not be appropriate for your analysis. Reading through this tutorial will give you a better understanding of how these (and other) metrics should be applied, so you can apply them effectively throughout your own explorations of these data.

Metrics covered in this tutorial:

• Firing rate (firing_rate)

• Presence ratio (presence_ratio)

• Amplitude cutoff (amplitude_cutoff)

• ISI violations (isi_violations)

• SNR (snr)

• Isolation distance (isolation_distance)

• d-prime (d_prime)

• Nearest-neighbors hit rate (nn_hit_rate)

How these metrics were calculated

The Python code used to calculate these metrics from the outputs of Kilosort2 is available in the ecephys_spike_sorting repository.

This code was incorporated into the SpikeMetrics repository by the SpikeInterface team. It’s now available as a PyPi package (pip install spikemetrics) if you’d like to try them out on your own data.

If you have any questions about the specific implementation of these metrics, or recommendations for new ones to include, we encourage you to submit an issue in either GitHub repository.

Accessing the metrics

Because these metrics are so important to interpreting your results, they are included in every DataFrame that stores information about individual units.

We will load the data from a Visual Coding Neuropixels session and look at the units table

import pynwb
import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt

#load a nwb file for an example ephys experiment
nwb_path = "/data/387858_2018-07-12_13-56-26_nwb_2026-08-19_07-52-55/387858_2018-07-12_13-56-26_nwb_2026-08-19_07-52-55.nwb.zarr"
nwbfile = pynwb.read_nwb(nwb_path)

units = nwbfile.units[:]
units.head()

	waveform_duration	cluster_id	peak_channel_id	cumulative_drift	amplitude_cutoff	snr	recovery_slope	isolation_distance	nn_miss_rate	silhouette_score	...	lifetime_sparseness_rf	lifetime_sparseness_sg	anterior_posterior_ccf_coordinate	dorsal_ventral_ccf_coordinate	left_right_ccf_coordinate	ecephys_structure_id	ecephys_structure_acronym	probe_horizontal_position	probe_vertical_position	ecephys_probe_id
id
950921187	0.604355	4	850249267	481.80	0.425574	2.209140	-0.118430	17.537571	0.009496	0.036369	...	0.068347	0.027059	7955.0	3756.0	7199.0	549.0	TH	11.0	20.0	729445650.0
950921172	0.521943	3	850249267	681.53	0.390098	1.959983	-0.109729	14.677643	0.003857	0.103446	...	0.226443	0.036345	7955.0	3756.0	7199.0	549.0	TH	11.0	20.0	729445650.0
950921152	0.467002	2	850249267	1070.71	0.500000	2.522905	-0.109867	15.783665	0.017776	0.027818	...	0.085703	0.030548	7955.0	3756.0	7199.0	549.0	TH	11.0	20.0	729445650.0
950921135	0.467002	1	850249267	253.42	0.500000	2.803475	-0.150379	26.666930	0.023742	0.076530	...	0.122004	0.011320	7955.0	3756.0	7199.0	549.0	TH	11.0	20.0	729445650.0
950921111	0.439531	0	850249267	141.82	0.018056	4.647943	-0.328727	66.901065	0.006595	NaN	...	0.030071	0.005186	7955.0	3756.0	7199.0	549.0	TH	11.0	20.0	729445650.0

5 rows × 94 columns

The full units table is unpacked more in dataset pages, but you can see that there are many columns, including several that pertain to quality control.

Note

The exact set of metrics in the units table will vary slightly between different datasets so look closely at the documentation for your dataset. At the same time, most of these quality metrics are used across our Neuropixels projects.

Metric name	Description
default_qc	True if the unit passes the default QC thresholds (see below)
isi_violations	rate of refractory-period violations, normalized (contamination measure; default threshold ≤ 0.5)
isi_violations_count	raw count of refractory-period violations
rp_violations	number of violations within a fixed refractory period
rp_contamination	estimated contamination from refractory-period violations (0–1)
sliding_rp_violation	minimum contamination estimated over a sliding refractory period, at ≥90% confidence (IBL method; 0–1)
snr	waveform signal-to-noise ratio on the peak channel
nn_hit_rate	nearest-neighbor hit rate; higher means better separated
sync_spike_2, sync_spike_4, sync_spike_8	fraction of spikes synchronous with 2 / 4 / 8 other units, flagging spikes shared across units
presence_ratio	fraction of the session with detected spikes (default threshold ≥ 0.8)
amplitude_cutoff	estimated fraction of spikes missed, from the amplitude distribution (default threshold ≤ 0.1)
nn_miss_rate	nearest-neighbor miss rate; estimates missing spikes
isolation_distance	distance to the nearest cluster in feature space (higher is better)
l_ratio	cluster-membership contamination measure (lower is better)
d_prime	linear-discriminability of this unit from its neighbors (higher is better)
silhouette	silhouette score of the unit’s cluster (higher is better)
amplitude_median	median spike amplitude (µV)
amplitude	spike amplitude (peak-to-trough)
amplitude_cv_median	median of the amplitude coefficient of variation over time
amplitude_cv_range	range of the amplitude coefficient of variation over time
drift_std	standard deviation of spike-position drift
drift_mad	median absolute deviation of drift
drift_ptp	peak-to-peak (maximum) drift

Three metrics contribute to the default_qc evaluation. These are isi_violation_ratio, presence_ratio, and amplitude_cutoff. The default filter values are as follows:

• isi_violations < 0.5

• amplitude_cutoff < 0.1

• presence_ratio > 0.9

Let’s look in more detail at the distribution of some quality metrics across the units in this Visual Coding Neuropixels asset. We’ll start by creating a function for plotting each metric in an aesthetically pleasing way:

from scipy.ndimage.filters import gaussian_filter1d
plt.rcParams.update({'font.size': 14})

def plot_metric(data, bins, x_axis_label, color, max_value=-1):

    h, b = np.histogram(data, bins=bins, density=True)

    x = b[:-1]
    y = gaussian_filter1d(h, 1)

    plt.plot(x, y, color=color)
    plt.xlabel(x_axis_label)
    plt.gca().get_yaxis().set_visible(False)
    [plt.gca().spines[loc].set_visible(False) for loc in ['right', 'top', 'left']]
    if max_value < np.max(y) * 1.1:
        max_value = np.max(y) * 1.1
    plt.ylim([0, max_value])

    return max_value

/tmp/ipykernel_104/773522037.py:1: DeprecationWarning: Please import `gaussian_filter1d` from the `scipy.ndimage` namespace; the `scipy.ndimage.filters` namespace is deprecated and will be removed in SciPy 2.0.0.
  from scipy.ndimage.filters import gaussian_filter1d

Firing rate

First, let’s take a look at firing rate, which is the most straightforward metric to compute. Firing rate is equal to the total number of spikes divided by the number of seconds in the recording. We’ll create a density plot of firing rate across all units in the dataset:

data = units['firing_rate']
bins = np.linspace(0,50,100)

max_value = plot_metric(data, bins, 'Firing rate (Hz)', 'red')

Since there are many units with low firing rates, let’s use a log scale instead:

data = np.log10(units['firing_rate'])
bins = np.linspace(-3,2,100)

max_value = plot_metric(data, bins, 'log$_{10}$ firing rate (Hz)', 'red')

Based on this plot, you can clearly see the approximately lognormal distribution of firing rates, which has been described previously [Buzsáki and Mizuseki, 2014]. However, there’s more weight on the lower tail of the distribution, which is likely due to some units missing spikes as a result of thresholding or drift. If we filter out contaminated units using another metric, nn_hit_rate (more on what this means later), the distribution becomes almost perfectly lognormal:

data = np.log10(units[units.nn_hit_rate > 0.9]['firing_rate'])
bins = np.linspace(-3,2,100)

max_value = plot_metric(data, bins, 'log$_{10}$ firing rate (Hz)', 'red')

Here’s a summary of things to keep in mind when using firing_rate in your analysis:

How it can be biased

• If a unit is poorly isolated, the firing rate will be over-estimated, because contaminating spikes will be included in the calculation

• If a unit’s amplitude is close to threshold, the firing rate will be under-estimated, because some spikes will be missing

• If a unit drifts out of the recording, the firing rate will be under-estimated, because spikes will not be detected for a portion of the recording

• If data acquisition is interrupted (true for a small subset of experiments), the firing rate will be under-estimated, because spikes will be missing from gaps in the recording

How it should be used

• Firing rate can be used to filter out units that have too few spikes to result in meaningful analysis. In this case, it may be better to use the firing rate for the specific interval you’re analyzing, because some units may drift out of the recording at other times.

• High firing rate units tend to be easier to isolate, since there are more spikes available for fitting the template in Kilosort2. However, there are other metrics that measure isolation more directly and would likely to be better to use instead.

Presence ratio

Presence ratio is not a standard metric in the field, but it’s straightforward to calculate and is an easy way to identify incomplete units. It measures the fraction of time during a session in which a unit is spiking, and ranges from 0 to 0.99 (an off-by-one error in the calculation ensures that it will never reach 1.0).

Let’s look at the distribution of presence ratio in this asset:

data = units['presence_ratio']
bins = np.linspace(0,1,100)
max_value = plot_metric(data, bins, 'Presence ratio', 'red')

It’s clear that most units have a presence ratio of 0.8 or higher, which means they are present for at least 80% of the recording. Units with lower presence ratio are likely to have drifted out of the recording, or had waveforms that changed so dramatically they were assigned to separate clusters.

Here’s a summary of things to keep in mind when using presence_ratio in your analysis:

How it can be biased

• Just because a unit has a high presence ratio doesn’t mean it’s immune to drift. If a unit’s amplitude drifts closer to the spike detection threshold, it can result in dramatic changes in apparent firing rate, even if the underlying physiology remains the same.

• Sometimes a low presence ratio can result from highly selective spiking patterns (e.g., firing only during running epochs)

How it should be used

• If you are analyzing changes in firing rate over the entire recording session, or are comparing responses to stimuli presented at the beginning and end of the experiment, presence ratio is a simple way to exclude units that would bias your results. However, you should also look at other quality metrics, such as amplitude cutoff, to check for more subtle effects of electrode drift.

• If you are only analyzing a short segment of the experiment, it may be helpful to disable the default presence ratio filter, in order to maximize the number of units available to you.

• If you’re unsure whether a unit has a low presence ratio due to electrode drift or selective firing, plotting its spike amplitudes over time can be informative.

Amplitude cutoff

Amplitude cutoff provides another way to check for units that are missing spikes. Unlike presence ratio, which detects units that drift out of the recording, amplitude cutoff provides an estimate of the false negative rate—e.g., the fraction of spikes below the spike detection threshold. Thus, amplitude cutoff is a measure of unit “completeness” that is complementary to presence ratio.

Let’s take a look at the distribution of values for amplitude cutoff in this asset:

data = units['amplitude_cutoff']
bins = np.linspace(0,0.5,200)
max_value = plot_metric(data, bins, 'Amplitude cutoff', 'red')

Amplitude cutoff is calculated from the distribution of spike amplitudes for each unit. This metric measures the degree to which this distribution is truncated, or “cut off,” as a proxy for the fraction of missing spikes. So an amplitude cutoff of, say, 0.1 would indicate that approximately 10% of spikes are missing from this unit.

If the peak of the amplitude distribution occurs at its lowest value, it’s impossible to estimate the fraction of missing spikes. In this case, the amplitude cutoff is set to 0.5. That explains why there are large peaks at both ends of the distribution, one around 0 and one at 0.5.

We can check the fraction of units with the maximum amplitude cutoff using the following code:

np.around(np.sum(units.amplitude_cutoff == 0.5) / len(units), 2)

0.21

Here’s a summary of things to keep in mind when using amplitude_cutoff in your analysis:

How it can be biased

• The calculation assumes that the amplitude histogram is symmetrical (i.e., it uses the upper tail of the distribution to estimate the fraction of spikes missing from the lower tail). If a unit’s waveform amplitude changes as a result of electrode drift, this assumption is usually invalid.

• Amplitude cutoff is only weakly correlated with other measures of unit quality, meaning it’s possible to have well-isolated units with high amplitude cutoff.

How it should be used

• If you are performing analyses that depends on precise measurements of spike timing, setting a low amplitude cutoff threshold (0.01 or lower) is recommended. This will remove a large fraction of units, but will ensure that the unit of interest contain most of the relevant spikes.

ISI violations

Inter-spike-interval (ISI) violations are a classic measure of unit contamination. Because all neurons have a biophysical refractory period, we can assume that any spikes occurring in rapid succession (<1.5 ms intervals) come from two different neurons. Therefore, the more a unit is contaminated by spikes from multiple neurons, the higher its isi_violations value will be.

The calculation for ISI violations comes from Hill et al. [2011]. Rather than reporting the fraction of spikes with ISI violations, their metric reports the relative firing rate of the hypothetical neurons that are generating these violations. You can interpret an ISI violations value of 0.5 as meaning that contaminating spikes are occurring at roughly half the rate of “true” spikes for that unit. In cases of highly contaminated units, the ISI violations value can sometimes be even greater than 1.

Let’s look at the distribution of ISI violations in this asset:

data = units['isi_violations']
bins = np.linspace(0,0.5,200)
max_value = plot_metric(data, bins, 'ISI violations', 'red')

Here’s a summary of things to keep in mind when using isi_violations in your analysis:

How it can be biased

• As with all metrics, ISI violations may not be stable throughout the experiment. It may be helpful to re-calculate it for the specific epochs you’re analyzing.

• Two neurons with similar waveforms, but firing in largely non-overlapping epochs, could end up being merged into the same cluster. In this case, the ISI violations may be low, even though the resulting unit is a highly contaminated. This situation would tricky to catch, but fortunately shouldn’t happen very often.

How it should be used

• Setting your ISI violations threshold to 0 (or close to it), will help ensure that contaminated units don’t make it into your analysis, but will greatly reduce the number of units available. You should think carefully about what degree of contamination your analysis can tolerate without biasing your conclusions. For example, if you are comparing firing rates of individual units across areas, you’ll want to set a low ISI violations threshold to prevent contaminating spikes from affecting your estimates. On the other hand, if you’re comparing overall firing rates between areas, counting spikes from contaminated clusters may be valid.

SNR

Signal-to-noise ratio, or SNR, is another classic metric of unit quality. It measures the ratio of the maximum amplitude of the mean spike waveform to the standard deviation of the background noise on one channel. Even though it’s widely used in the literature, we don’t recommend using it on Neuropixels data for two reasons:

1. It only takes into account the unit’s peak channel, despite the fact that waveforms are often spread across a dozen channels or more.

2. If the waveform changes due to drift, peak channel SNR can change dramatically, even though overall isolation quality remains consistent.

Nevertheless, it can still be helpful to look at the distribution of SNRs:

data = units['snr']
bins = np.linspace(0,10,100)
max_value = plot_metric(data, bins, 'SNR', 'red')

Here’s a summary of things to keep in mind when using snr in your analysis:

How it can be biased

• SNR only considers information contained in a single channel, and therefore cannot capture isolation quality accurately

• If the peak channel moves across the probe due to electrode drift, SNR will drop.

How it should be used

• SNR can be helpful as a point of comparison to the previous literature

• It should not be used to filter data in isolation

That said, a modified version of the SNR metric that is tolerant to electrode drift could be highly informative. Future Allen Institute data releases may include such a metric.

Isolation distance

Isolation distance is a metric based on the principal components (PCs) of a unit’s waveforms. After the spike sorting step is complete, the waveforms for every spike are projected into a lower-dimensional principal component space. By default, Kilosort2 saves the top 3 PCs for 32 channels around each unit’s peak channel—this is a huge amount of data, but it’s greatly compressed compared to the original 60 samples x 350+ channels for each waveform. PC-based metrics are a useful way of validating cluster quality because, at least for Kilosort2, the original sorting process doesn’t rely on the waveform’s principal components.

You can imagine each unit’s PCs a clusters in a 32 x 3 = 96-dimensional space. Isolation distance calculates the size of the 96-dimensional sphere that includes as many “other” spikes as are contained in the original unit’s cluster, after normalizing the clusters by their standard deviation in each dimension (Mahalanobis distance). The higher the isolation distance, the more a unit is separated from its neighbors in PC space, and therefore the lower the likelihood that it’s contaminated by spikes from multiple units.

Let’s look at the range of isolation distances in this asset:

bins = np.linspace(0,170,50)
data = units['isolation_distance']
max_value = plot_metric(data, bins, 'Isolation distance', 'red')

Here’s a summary of things to keep in mind when using isolation_distance in your analysis:

How it can be biased

• Isolation distance is not immune to drift; if a unit’s waveform changes as a result of electrode motion, it could reduce isolation distance without necessarily causing the unit to become more contaminated.

• The exact value of isolation distance will depend on the number of PCs used in the calculation; therefore, it’s difficult to compare this metric to previous reports in the literature.

How it should be used

• Isolation distance is correlated with overall cluster quality, but it’s not a direct measure of contamination rate. For this reason, it should be used in conjunction with other metrics, such as isi_violations, that more directly measure the likelihood of contaminating spikes.

d-prime

Like isolation distance, d-prime is another metric calculated for the waveform PCs. It uses linear discriminant analysis to calculate the separability of one unit’s PC cluster and all of the others. A higher d-prime value indicates that the unit is better isolated from its neighbors.

bins = np.linspace(0,15,50)
data = units['d_prime']
max_value = plot_metric(data, bins, 'd-prime', 'red')

Here’s a summary of things to keep in mind when using d_prime in your analysis:

How it can be biased

• Like isolation distance, d-prime is not tolerant to drift. Since a single value of d-prime is computed for the entire session, the d-prime value is actually a lower bound on the true value of this metric computed at any one timepoint.

How it should be used

• d-prime, in principal, gives you an estimate of the false positive rate for each unit. However, more work is required to validate this.

Nearest-neighbors hit rate

Nearest-neighbors hit rate is another PC-based quality metric. It’s derived from the ‘isolation’ metric originally reported in Chung et al. [2017]. This metric looks at the PCs for one unit and calculates the fraction of their nearest neighbors that fall within the same cluster. If a unit is highly contaminated, then many of the closest spikes will come from other units. Nearest-neighbors hit rate is nice because it always falls between 0 and 1, making it straightforward to compare across different datasets.

bins = np.linspace(0,1,100)
data = units['nn_hit_rate']
max_value = plot_metric(data, bins, 'Nearest-neighbors hit rate', 'red')

Here’s a summary of things to keep in mind when using nn_hit_rate in your analysis:

How it can be biased

• Like the other PC-based metrics, nn_hit_rate can be negatively impacted by electrode drift.

How it should be used

• nn_hit_rate is a nice proxy for overall cluster quality, but should be used in conjunction with other metrics that measure missing spikes or contamination rate more directly.

Summary

To summarize, let’s take a look at the range of values that each of these metrics takes across the whole Visual Coding dataset:

metrics = ['firing_rate',
           'presence_ratio',
           'amplitude_cutoff',
           'isi_violations',
           'snr',
           'isolation_distance',
           'd_prime',
           'nn_hit_rate']

ranges = [[0,20],
          [0.9,0.995],
          [0,0.5],
          [0,2],
          [0,8],
          [0,125],
          [0,10],
          [0,1]]

_ = plt.figure(figsize=(5,10))

for idx, metric in enumerate(metrics):

    data = units[metric].values
    data = data[np.invert(np.isnan(data))]

    _ = plt.subplot(len(metrics),1,idx+1)
    _ = plt.boxplot(data, showfliers=False, showcaps=False, vert=False)
    _ = plt.ylim([0.8,1.2])
    _ = plt.xlim(ranges[idx])
    _ = plt.yticks([])

    plt.title(metric)

plt.tight_layout()