Lussac modules

Here, we describe each of Lussac modules, what they do, how they do it, and what parameters are available.

How waveforms are extracted

In many different modules, waveforms need to be extracted a wvf_extraction parameter dictionary is used to know how waveforms should be extracted and processed. The different parameters are:

ms_before: how many samples (in ms) should be extracted before the spike peak.
ms_after: how many samples (in ms) should be extracted after the spike peak.
max_spikes_per_unit: how many spikes should be extracted per unit to create the template.
filter_band: a list of 2 floats, representing the minimum and maximum frequency (in Hz) for a Gaussian bandpass filter. Note that Lussac extracts an extra margin to remove border filtering artifacts.

The `units_categorization` module

This module will label units as belonging to a certain category if they meet some criteria. If a unit already belongs to a category, it will not be re-categorized. So the order matters!

This module takes as a key the name of the category, and as a value a dictionary containing the criteria. Each criterion return a value for each unit, and a minimum and/or maximum can be set.

You can also specify the parameters for wvf_extraction (max_spikes_per_unit, ms_before, ms_after, filter).

firing_rate: returns the mean firing rate of the unit (in Hz).
contamination: returns the estimated contamination of the unit (between 0 and 1; 0 being pure). The refractory_period = [censored_period, refractory_period] has to be set (in ms).
amplitude: returns the mean amplitude of the unit’s template (in µV if the recording object can be scaled to µV). Optional parameters can be set to the function spikeinterface.core.get_template_extremum_amplitude.
SNR: returns the signal-to-noise ratio of the unit. Optional parameters can be set to the function spikeinterface.qualitymetrics.compute_snrs
sd_ratio: returns the standard deviation in the spike amplitudes for the unit divided by the standard deviation on the same channel. Optional parameters can be set to the function spikeinterface.postprocessing.compute_spike_amplitudes as a dictionary with key spike_amplitudes_kwargs. Optional parameters can be set to the function spikeinterface.qualitymetrics.compute_sd_ratio as a dictionary with key sd_ratio_kwargs.
ISI_portion: Returns the fraction (between 0 and 1) of inter-spike intervals inside a time range. the range = [min_t, max_t] has to be set (in ms).

Example of categorization

Here is an example for categorizing complex-spikes from more regular spikes (cerebellar cortex example):

"units_categorization": {
    "all": {  // Categorize all units.
        "wvf_extraction": {  // Parameters for the waveform extraction.
            "ms_before": 1.0,
            "ms_after": 1.5,
            "max_spikes_per_unit": 500,
            "filter": [150.0, 7000.0]  // Gaussian bandpass filter with cutoffs at 150 and 7,000 Hz.
        },
        "CS": {  // Criteria for complex-spikes category.
            "firing_rate": {  // Firing rate < 5 Hz
                "max": 5.0
            },
            "ISI_portion": {  // Few spikes between 10 and 35 ms in the ISI.
                "range": [10.0, 35.0],
                "max": 0.05
            }
        },
        "spikes": {  // Categorize more "regular" spikes
            "firing_rate": {
                "min": 0.4,
                "max": 200.0
            },
            "contamination": {  // Maximum of 30% contamination.
                "refractory_period": [0.3, 1.0],
                "max": 0.3
            },
            "SNR": {
                "peak_sign": "neg",  // Example of an SI parameter.
                "min": 2.5
            }
        }
    }
}

Clearing units category

It is possible to remove the category label on units, by setting the category name to "clear". For example:

"units_categorization": {
    "all": {"clear": {}}  // Clear category label from all units.
}

The `align_units` module

This module will align units by using their template waveform. The algorithm is not that straightforward:

First, a threshold is set and we look at the first peak that is higher than this threshold (both in positive and negative values). Next, the algorithm checks a few more samples in time for a higher peak (if one is detected, it takes this one).

The rational behind it is that the “center” of the spike should be when the neuron starts its action potential. For multi-phasic spikes, this is usually the first one. The threshold and check_next parameters are here to make sure we’re not detecting noise.

TODO: Insert image with example of sub-threshold peak and check_next.

This module’s parameters are:

wvf_extraction: to construct the templates. The ms_before and ms_after parameters determine the max shift for alignment.
threshold (optional): Threshold multiplicator (between 0 and 1). The real threshold is max(template) * threshold. By default: 0.5
check_next (optional): Number of samples to check after the first peak (put 0 to not check after the first peak). By default: 10

Example of units alignment

"align_units": {
    "all": {  // Align all units.
        "wvf_extraction": {
            "ms_before": 1.0,
            "ms_after": 2.0,
            "max_spikes_per_unit": 2000,  // Use 2,000 random spikes to construct templates.
            "filter_band": [300.0, 6000.0]  // Gaussian-filter between 300 and 6000 Hz.
        },
        "threshold": 0.5,  // Threshold at 50% of the maximum.
        "check_next": 5  // Check the next 5 samples.
    }
}

The `remove_bad_units` module

This module will remove the units that meet at least one of the criteria. The criteria are the same as those described in units_categorization.

Example of units removal

"remove_bad_units": {
    "CS": {  // Remove complex-spike units with contamination > 35%
        "wvf_extraction": {},  // If you want to change how the waveforms are extracted.
        "contamination": {
            "refractory_period": [1.5, 25.0],
            "max": 0.35
        }
    },
    "spikes": {  // Remove units with firing rate < 1.0 Hz or amplitude std > 80 µV
        "firing_rate": {
            "min": 1.0
        },
        "sd_ratio": {
            "max": 2.0
        }
    }

The `remove_duplicated_spikes` module

This module will remove spikes that are considered duplicates (i.e. too close to one another).

This is done by setting a censored_period window under which there cannot be 2 spikes.
Be careful! This is different from the refractory_period! It’s very useful to keep spikes in the refractory period to estimate the contamination. The censored period is designed to remove duplicated spikes.
Typical values of censored_period usually lie between 0.2 and 0.4 ms, whereas the refractory period is almost always greater than 0.9ms.

This module’s parameters are:

censored_period: in ms (by default, 0.3).
method (optional): method used to remove duplicates (used by spikeinterface.curation.find_duplicated_spikes). By default: "keep_first_iterative"

Example of duplicated spikes removal

"remove_duplicated_spikes": {
    "all": {
        "censored_period": 0.3
    }
}

The `remove_redundant_units` module

This module will look for redundant units in analyses (by looking at the rate of coincident spikes between units in individual analyses).

If redundant units are detected, all but one will be removed (the chosen one depends on the remove_strategy used).

This module’s parameters are:

wvf_extraction: to construct the templates (required depending on the remove strategy). If not required, just set it to null.
arguments: a dict containing the parameters to give to spikeinterface.curation.remove_redundant_units.

Example of redundant units removal

"remove_redundant_units": {
    "all": {
        "wvf_extraction": {
            "ms_before": 1.0,
            "ms_after": 1.5,
            "max_spikes_per_unit": 500
        },
        "arguments": {
            "align": true,  // Can be set to 'false' if you already used the 'align_units' module.
            "delta_time": 0.3,  // Window (in ms) to consider coincident spikes.
            "duplicate_threshold": 0.7,  // If coincidence >= 70%, consider the units redundant.
            "remove_strategy": "highest_amplitude"  // Keep the unit with the highest amplitude.
        }
    }
}

The `merge_units` module

This module looks for units that correspond to the same neuron (inside each individual analysis separately), and merges them together if the merge is deemed beneficial.

This is done by first looking over all pairs of units, and estimating if they likely come from the same neuron, on the basis of: proximity, matching correlograms, matching templates.
Then, pairs that don’t increase the quality score if the merge is performed are discarded. With this discard, the worse of both units is removed (because it usually is a bad split unit).
Finally, a graph is constructed from the remaining pairs. For each connected component (i.e. each putative neuron), we iteratively merge the best pair until everything is merged or there are no more merges that increase the quality score metric. If some unmerged units remain, they are discarded.

This modules parameters are:

refractory_period = [censored_period, refractory_period]: in ms. By default: [0.2, 1.0].
wvf_extraction: to construct the templates.
correlograms: a dict containing the parameters to construct the correlograms.
- window_ms: The total window size of the correlogram (in ms). A value of 100.0 will create a correlogram of size [-50.0, 50.0] ms. By default: 150 ms.
- bin_ms: The size of the bins in the correlogram (in ms). By default: 0.04 ms.
auto_merge_params: a dict containing the parameters to give to spikeinterface.curation.auto_merge.compute_merge_unit_groups.

Example of merging units

"merge_units": {
    "all": {
        "refractory_period": [0.2, 1.0],
        "wvf_extraction": {
            "ms_before": 1.0,
            "ms_after": 1.5,
            "max_spikes_per_unit": 2000,
            "filter_band": [150, 7000]
        },
        "correlograms": {
            "window_ms": 150,
            "bin_ms": 0.04
        },
        "auto_merge_params": {
            "steps_params": {
                "correlogram": {
                    "corr_diff_thresh": 0.16,
                    "censor_correlograms_ms": 0.2,
                    "sigma_smooth_ms": 0.6,
                    "adaptive_window_thresh": 0.5
                },
                "template_similarity": {"template_diff_thresh": 0.25}
            },
            "firing_contamination_balance": 2.5,  // k = 2.5 in the paper.
            "resolve_graph": false  // False by default because Lussac implements its own graph resolution.
        }
    }
}

The `merge_sortings` module

This module is the heart of Lussac. It will merge all individual analyses into a single one, following a complex algorithm.

STEP 1: Create a graph where each node is a unit, and each edge links similar units (based on the correlation of their spike trains).
STEP 2: Detect and remove merged units.
STEP 3: Detect “wrong” edges and remove them.
STEP 4: For each community, create/select the best unit.

The parameters used in this module are:

refractory_period = [censored_period, refractory_period]: in ms. By default: [0.2, 1.0].
max_shift: The maximum time shift when re-aligning pairs of units (in ms). By default: 1.33 ms.
require_multiple_sortings_match: Whether to remove lone units (i.e. units that are not matched with any other unit). By default: True.
similarity: a dict to compute the similarity (i.e. spike trains correlation) in STEP 1.
- min_similarity: The minimum similarity to consider two units similar. By default: 0.3.
- window: The maximum lag (in ms) allowed between two spikes to be considered similar. By default: 0.2 ms.
correlogram_validation: a dict to compute the validation correlogram in STEP 3.
- max_time: The maximum time for the correlogram (in ms). By default: 70 ms (i.e. correlogram computed between [-70, 70] ms).
- gaussian_std: The standard deviation of the Gaussian kernel used to smooth the correlogram (in ms). By default: 0.6 ms.
- gaussian_truncate: The Gaussian is truncated after X standard deviations for faster computation. By default: X = 5.
- bin_ms (optional): The size of the bins in the correlogram (in ms). By default: very small, adaptative to max_time.
waveform_validation: a dict to compute the validation waveform in STEP 3.
- wvf_extraction: to construct the templates. By default ms_before = 1.0, ms_after = 2.0, max_spikes_per_unit = 1000, filter_band = [250, 6000].
- num_channels: The number of channels used to compare waveforms. By default: 5.
merge_check: a dict to compute the merge check in STEP 2.
- cross_cont_threshold: The threshold above which the cross-contamination is considered too high. By default: 0.10. Note that the cross-contamination needs to be significantly higher, using a statistical test.
clean_edges: a dict with the thresholds used for STEP 3.
- template_diff_threshold: The threshold above which the template difference is considered too high. By default: 0.10.
- corr_diff_threshold: The threshold above which the correlation difference is considered too high. By default: 0.12.
- cross_cont_threshold: The threshold above which the cross-contamination is considered too high. By default: 0.06. Note that the cross-contamination needs to be significantly higher, using a statistical test.

Example of merging sortings

"merge_sortings": {
    "all": {
        "refractory_period": [0.2, 1.0],
        "max_shift": 1.33,
        "require_multiple_sortings_match": true,
        "similarity": {
            "min_similarity": 0.3,
            "window": 0.2
        },
        "correlogram_validation": {
            "max_time": 70.0,
            "gaussian_std": 0.6,
            "gaussian_truncate": 5.0
        },
        "waveform_validation": {
            "wvf_extraction": {
                "ms_before": 1.0,
                "ms_after": 2.0,
                "max_spikes_per_unit": 1000,
                "filter_band": [250.0, 6000.0]
            },
            "num_channels": 5
        },
        "merge_check": {
            "cross_cont_threshold": 0.10
        },
        "clean_edges": {
            "template_diff_threshold": 0.10,
            "corr_diff_threshold": 0.12,
            "cross_cont_threshold": 0.06
        }
    }
}

The `find_purkinje_cells` module

This module is only meant for cerebellar cortex recordings. It will link simple spikes and complex spikes coming from the same Purkinje cell, and set it as a property lussac_purkinje (this property is automatically exported in the export_to_phy module).

TODO: Explain how it works.

This module’s parameters are:

cross_corr_pause: the band over which to look for the pause (in ms). By default: [0.0, 8.0].
threshold: TODO
ss_min_fr: Minimum firing rate to consider putative simple spikes (in Hz). By default: 40.0.
cs_min_fr: Minimum firing rate to consider putative complex spikes (in Hz). By default: 0.5.
cs_max_fr: Maximum firing rate to consider putative complex spikes (in Hz). By default: 3.0.

Example of finding Purkinje cells

"find_purkinje_cells": {
    "all": {
        "cross_corr_pause": [0.0, 8.0],
        "threshold": 0.4,
        "ss_min_fr": 40.0,
        "cs_max_fr": 3.0
    }
}

The `export_to_phy` module

This module will export all sortings in their current state to the phy format (if merge_sortings was called before, will only export the merged sorting).

This module’s parameters are:

path: path to the folder where to export the sorting(s). If multiple sortings exists, a subfolder will be created for each of them.
wvf_extraction: to construct the templates.
export_params: a dict containing the parameters to give to spikeinterface.exporters.export_to_phy.
estimate_contamination (optional): a dict containing the refractory period for each category. If given, will output the estimated contamination of the units.

Example of export to phy

"export_to_phy": {
    "all": {
        "path": "$PARAMS_FOLDER/lussac/final_output",
        "wvf_extraction": {
            "ms_before": 1.0,
            "ms_after": 3.0,
            "max_spikes_per_unit": 1000
        },
        "export_params": {
            "compute_amplitudes": true,
            "compute_pc_features": false,
            "copy_binary": false,
            "template_mode": "average",
            "sparsity": {
                "method": "radius",
                "radius_um": 75.0
            },
            "verbose": false
        },
        "estimate_contamination": {
            "all": [0.3, 1.0]
        }
    }
}

The `export_to_sigui` module

This module will export all sortings in their current state to the SpikeInterface GUI format (if merge_sortings was called before, will only export the merged sorting).

This is equivalent to just a SortingAnalyzer with some extra arguments.

This module’s parameters are:

path: path to the folder where to export the sorting(s). If multiple sortings exists, a subfolder will be created for each of them.
wvf_extraction: to construct the templates.
spike_amplitudes (optional): either a dict or False. If a dict, will compute and export the spike amplitudes, the content of the dictionary being the parameters for spikeinterface.postprocessing.compute_spike_amplitudes. By default dict().
principal_components (optional): either a dict or False. If a dict, will compute and export the PCA, the content of the dictionary being the parameters for spikeinterface.postprocessing.compute_principal_components. By default False.

Example of export to SI GUI

"export_to_sigui": {
    "all": {
        "path": "$PARAMS_FOLDER/lussac/final_output",
        "wvf_extraction": {
            "ms_before": 1.0,
            "ms_after": 3.0,
            "max_spikes_per_unit": 1000
        }
    }
}

Lussac modules

How waveforms are extracted

The units_categorization module

Example of categorization

Clearing units category

The align_units module

Example of units alignment

The remove_bad_units module

Example of units removal

The remove_duplicated_spikes module

Example of duplicated spikes removal

The remove_redundant_units module

Example of redundant units removal

The merge_units module

Example of merging units

The merge_sortings module

Example of merging sortings

The find_purkinje_cells module

Example of finding Purkinje cells

The export_to_phy module

Example of export to phy

The export_to_sigui module

Example of export to SI GUI

The `units_categorization` module

The `align_units` module

The `remove_bad_units` module

The `remove_duplicated_spikes` module

The `remove_redundant_units` module

The `merge_units` module

The `merge_sortings` module

The `find_purkinje_cells` module

The `export_to_phy` module

The `export_to_sigui` module