Neural circuits encode prior knowledge of temporal statistics.

1. from typing import Any
2. import numpy as np
3. from overrides import override
4. from lussac.core import MonoSortingModule
5. import spikeinterface.core as si
6. import spikeinterface.postprocessing as spost
7. import spikeinterface.metrics as sm
10. class FindPurkinjeCells(MonoSortingModule):
11. """
12. This module is meant purely for the cerebellar cortex.
13. Purkinje cells are big GABAergic neurons that can fire both simple spikes and complex spikes.
14. Since simple spikes and complex spikes have different shape, they will be clustered in different units.
15. This module will find simple spikes and complex spikes coming from the same Purkinje cell,
16. and label then into a 'lussac_purkinje' property in the sorting object.
17. This is done by looking at the cross-correlogram between the two units: after a complex spike,
18. there is always a pause of at least 8-10 ms before simple spike activity resumes.
19. Also, when high-pass filtered aggressively, the templates are nearly identical.
20. """
22. @property
23. @override
24. def default_params(self) -> dict[str, Any]:
25. return {
26. 'cross_corr_pause': [0.0, 8.0],
27. 'threshold': 0.4,
28. 'ss_min_fr': 40.0,
29. 'cs_min_fr': 0.5,
30. 'cs_max_fr': 3.0
31. }
33. @override
34. def update_params(self, params: dict[str, Any]) -> dict[str, Any]:
35. params = super().update_params(params)
36. if params['ss_min_fr'] <= params['cs_max_fr']:
37. raise ValueError("Error in 'find_purkinje_cells': 'ss_min_fr' must be greater than 'cs_max_fr'")
38. if params['threshold'] < 0.0:
39. raise ValueError("Error in 'find_purkinje_cells': 'threshold' cannot be negative")
41. return params
43. @override
44. def run(self, params: dict[str, Any]) -> si.BaseSorting:
45. analyzer = si.SortingAnalyzer.create(self.sorting, self.recording)
46. firing_rates = sm.compute_firing_rates(analyzer)
48. putative_ss_units = [unit_id for unit_id, mean_fr in firing_rates.items() if mean_fr >= params['ss_min_fr']]
49. putative_cs_units = [unit_id for unit_id, mean_fr in firing_rates.items() if params['cs_min_fr'] <= mean_fr <= params['cs_max_fr']]
51. sorting = self.sorting.select_units([*putative_ss_units, *putative_cs_units])
52. correlograms, bins = spost.compute_correlograms(sorting, window_ms=25.0, bin_ms=1.0, method="numba")
53. mask = ((bins >= params['cross_corr_pause'][0]) & (bins < params['cross_corr_pause'][1]))[:-1]
55. ss_cs_pairs = []
56. for ss_id in putative_ss_units:
57. for cs_id in putative_cs_units:
58. ss_ind = sorting.id_to_index(ss_id)
59. cs_ind = sorting.id_to_index(cs_id)
60. cross_corr = correlograms[ss_ind, cs_ind, :]
62. baseline = np.median(cross_corr[bins[:-1] < 0.0])
63. if np.median(cross_corr[mask]) < baseline * params['threshold']: # Check for pause.
64. if np.median(cross_corr[mask[::-1]] < baseline * params['threshold']): # Check for asymmetry.
65. continue
66. ss_cs_pairs.append((ss_id, cs_id)) # TODO: Also check templates?
68. lussac_purkinje = np.empty(self.sorting.get_num_units(), dtype=object)
69. paired_units = np.unique(np.array(ss_cs_pairs).flatten())
70. for unit_id in paired_units:
71. unit_ind = self.sorting.id_to_index(unit_id)
72. mask = np.array([unit_id in ss_cs_pairs[i] for i in range(len(ss_cs_pairs))])
73. text = " ; ".join([f"{ss_id}-{cs_id}" for ss_id, cs_id in np.array(ss_cs_pairs)[mask]])
74. lussac_purkinje[unit_ind] = text
76. self.sorting.set_property('lussac_purkinje', lussac_purkinje)
77. return self.sorting

find_purkinje_cells.py at commit 2a0dce1, under AGPL-3.0 · at the source