OSCR

Sequential visual stimuli increase high frequency power in the visual cortex.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § Materials and methods › Multi-unit and Single-unit extraction and analysis ↔ ecephys_spike_sorting/modules/quality_metrics/metrics.py, lines 20–157 · score 0.80 · quality metrics, isi violation, Isolation distance, firing rate, pre, clustering

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 997 lines · 35 KB · BSD-2-Clause · 1 match

  1. import numpy as np
  2. import pandas as pd
  3. from collections import OrderedDict
  5. import warnings
  7. from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
  8. from sklearn.neighbors import NearestNeighbors
  9. from sklearn.metrics import silhouette_score
  11. from scipy.spatial.distance import cdist
  12. from scipy.stats import chi2
  13. from scipy.ndimage.filters import gaussian_filter1d
  15. from ...common.epoch import Epoch
  16. from ...common.utils import printProgressBar, get_spike_depths
  20. def calculate_metrics(spike_times,
  21. spike_clusters,
  22. spike_templates,
  23. amplitudes,
  24. channel_map,
  25. pc_features,
  26. pc_feature_ind,
  27. params,
  28. epochs = None):
  30. """ Calculate metrics for all units on one probe
  32. Inputs:
  33. ------
  34. spike_times : numpy.ndarray (num_spikes x 0)
  35. Spike times in seconds (same timebase as epochs)
  36. spike_clusters : numpy.ndarray (num_spikes x 0)
  37. Cluster IDs for each spike
  38. spike_templates : numpy.ndarray (num_spikes x 0)
  39. Original template IDs for each spike time
  40. amplitudes : numpy.ndarray (num_spikes x 0)
  41. Amplitude value for each spike time
  42. channel_map : numpy.ndarray (num_units x 0)
  43. Original data channel for pc_feature_ind array
  44. pc_features : numpy.ndarray (num_spikes x num_pcs x num_channels)
  45. Pre-computed PCs for blocks of channels around each spike
  46. If 'None', PC-based metrics will not be computed
  47. pc_feature_ind : numpy.ndarray (num_units x num_channels)
  48. Channel indices of PCs for each unit
  49. epochs : list of Epoch objects
  50. contains information on Epoch start and stop times
  51. params : dict of parameters
  52. 'isi_threshold' : minimum time for isi violations
  55. Outputs:
  56. --------
  57. metrics : pandas.DataFrame
  58. one column for each metric
  59. one row per unit per epoch
  61. """
  63. metrics = pd.DataFrame()
  65. np.random.seed(9999)
  67. if epochs is None:
  68. epochs = [Epoch('complete_session', 0, np.inf)]
  70. total_units = len(np.unique(spike_clusters))
  71. total_epochs = len(epochs)
  73. for epoch in epochs:
  75. in_epoch = (spike_times > epoch.start_time) * (spike_times < epoch.end_time)
  77. print("Calculating isi violations")
  78. isi_viol = calculate_isi_violations(spike_times[in_epoch], spike_clusters[in_epoch], total_units, params['isi_threshold'], params['min_isi'])
  80. print("Calculating corrected isi violations")
  81. isi_viol_corrected = calculate_isi_violations_corrected(spike_times[in_epoch], spike_clusters[in_epoch], total_units, params['isi_threshold'], params['min_isi'])
  83. print("Calculating presence ratio")
  84. presence_ratio = calculate_presence_ratio(spike_times[in_epoch], spike_clusters[in_epoch], total_units)
  86. print("Calculating firing rate")
  87. firing_rate = calculate_firing_rate(spike_times[in_epoch], spike_clusters[in_epoch], total_units)
  89. print("Calculating amplitude cutoff")
  90. amplitude_cutoff = calculate_amplitude_cutoff(spike_clusters[in_epoch], amplitudes[in_epoch], total_units)
  92. if pc_features is not None:
  94. print("Calculating PC-based metrics")
  95. isolation_distance, l_ratio, d_prime, nn_hit_rate, nn_miss_rate = calculate_pc_metrics(spike_clusters[in_epoch],
  96. spike_templates[in_epoch],
  97. total_units,
  98. pc_features[in_epoch,:,:],
  99. pc_feature_ind,
  100. params['num_channels_to_compare'],
  101. params['max_spikes_for_unit'],
  102. params['max_spikes_for_nn'],
  103. params['n_neighbors'])
  107. print("Calculating silhouette score")
  108. the_silhouette_score = calculate_silhouette_score(spike_clusters[in_epoch],
  109. spike_templates[in_epoch],
  110. total_units,
  111. pc_features[in_epoch,:,:],
  112. pc_feature_ind,
  113. params['n_silhouette'])
  115. print("Calculating drift metrics")
  116. max_drift, cumulative_drift = calculate_drift_metrics(spike_times[in_epoch],
  117. spike_clusters[in_epoch],
  118. spike_templates[in_epoch],
  119. total_units,
  120. pc_features[in_epoch,:,:],
  121. pc_feature_ind,
  122. params['drift_metrics_interval_s'],
  123. params['drift_metrics_min_spikes_per_interval'])
  125. cluster_ids = np.unique(spike_clusters)
  127. epoch_name = [epoch.name] * len(cluster_ids)
  129. if pc_features is not None:
  130. metrics = pd.concat((metrics, pd.DataFrame(data= OrderedDict((('cluster_id', cluster_ids),
  131. ('firing_rate' , firing_rate),
  132. ('presence_ratio' , presence_ratio),
  133. ('isi_viol' , isi_viol),
  134. ('isi_viol_corrected' , isi_viol_corrected),
  135. ('amplitude_cutoff' , amplitude_cutoff),
  136. ('isolation_distance' , isolation_distance),
  137. ('l_ratio' , l_ratio),
  138. ('d_prime' , d_prime),
  139. ('nn_hit_rate' , nn_hit_rate),
  140. ('nn_miss_rate' , nn_miss_rate),
  141. ('silhouette_score', the_silhouette_score),
  142. ('max_drift', max_drift),
  143. ('cumulative_drift', cumulative_drift),
  144. ('epoch_name' , epoch_name),
  145. )))))
  147. else:
  148. metrics = pd.concat((metrics, pd.DataFrame(data= OrderedDict((('cluster_id', cluster_ids),
  149. ('firing_rate' , firing_rate),
  150. ('presence_ratio' , presence_ratio),
  151. ('isi_viol' , isi_viol),
  152. ('isi_viol_corrected' , isi_viol_corrected),
  153. ('amplitude_cutoff' , amplitude_cutoff),
  154. ('epoch_name' , epoch_name),
  155. )))))
  157. return metrics
  159. # ===============================================================
  161. # HELPER FUNCTIONS TO LOOP THROUGH CLUSTERS:
  163. # ===============================================================
  165. def calculate_isi_violations(spike_times, spike_clusters, total_units, isi_threshold, min_isi):
  167. cluster_ids = np.unique(spike_clusters)
  169. viol_rates = np.zeros((total_units,))
  171. for idx, cluster_id in enumerate(cluster_ids):
  173. printProgressBar(idx + 1, total_units)
  175. for_this_cluster = (spike_clusters == cluster_id)
  176. viol_rates[idx], num_violations = isi_violations(spike_times[for_this_cluster],
  177. min_time = np.min(spike_times),
  178. max_time = np.max(spike_times),
  179. isi_threshold=isi_threshold,
  180. min_isi = min_isi)
  182. return viol_rates
  184. def calculate_isi_violations_corrected(spike_times, spike_clusters, total_units, isi_threshold, min_isi):
  186. cluster_ids = np.unique(spike_clusters)
  188. viol_rates = np.zeros((total_units,))
  190. for idx, cluster_id in enumerate(cluster_ids):
  192. printProgressBar(idx + 1, total_units)
  194. for_this_cluster = (spike_clusters == cluster_id)
  195. viol_rates[idx], num_violations = isi_violations_corrected(spike_times[for_this_cluster],
  196. min_time = np.min(spike_times),
  197. max_time = np.max(spike_times),
  198. isi_threshold=isi_threshold,
  199. min_isi = min_isi)
  201. return viol_rates
  203. def calculate_presence_ratio(spike_times, spike_clusters, total_units):
  205. cluster_ids = np.unique(spike_clusters)
  207. ratios = np.zeros((total_units,))
  209. for idx, cluster_id in enumerate(cluster_ids):
  211. printProgressBar(idx + 1, total_units)
  213. for_this_cluster = (spike_clusters == cluster_id)
  214. ratios[idx] = presence_ratio(spike_times[for_this_cluster],
  215. min_time = np.min(spike_times),
  216. max_time = np.max(spike_times))
  218. return ratios
  222. def calculate_firing_rate(spike_times, spike_clusters, total_units):
  224. cluster_ids = np.unique(spike_clusters)
  226. firing_rates = np.zeros((total_units,))
  228. min_time = np.min(spike_times)
  229. max_time = np.max(spike_times)
  231. for idx, cluster_id in enumerate(cluster_ids):
  233. printProgressBar(idx + 1, total_units)
  235. for_this_cluster = (spike_clusters == cluster_id)
  236. firing_rates[idx] = firing_rate(spike_times[for_this_cluster],
  237. min_time = np.min(spike_times),
  238. max_time = np.max(spike_times))
  240. return firing_rates
  243. def calculate_amplitude_cutoff(spike_clusters, amplitudes, total_units):
  245. cluster_ids = np.unique(spike_clusters)
  247. amplitude_cutoffs = np.zeros((total_units,))
  249. for idx, cluster_id in enumerate(cluster_ids):
  251. printProgressBar(idx + 1, total_units)
  254. for_this_cluster = (spike_clusters == cluster_id)
  255. amplitude_cutoffs[idx] = amplitude_cutoff(amplitudes[for_this_cluster])
  257. return amplitude_cutoffs
  260. def calculate_pc_metrics_one_cluster(cluster_peak_channels, idx, cluster_id,cluster_ids,
  261. half_spread, pc_features, pc_feature_ind,
  262. spike_clusters, spike_templates,
  263. max_spikes_for_cluster, max_spikes_for_nn, n_neighbors):
  265. peak_channel = cluster_peak_channels[idx]
  266. num_spikes_in_cluster = np.sum(spike_clusters == cluster_id)
  268. half_spread_down = peak_channel \
  269. if peak_channel < half_spread \
  270. else half_spread
  272. half_spread_up = np.max(pc_feature_ind) - peak_channel \
  273. if peak_channel + half_spread > np.max(pc_feature_ind) \
  274. else half_spread
  276. channels_to_use = np.arange(peak_channel - half_spread_down, peak_channel + half_spread_up + 1)
  277. units_in_range = cluster_ids[np.isin(cluster_peak_channels, channels_to_use)]
  279. spike_counts = np.zeros(units_in_range.shape)
  281. for idx2, cluster_id2 in enumerate(units_in_range):
  282. spike_counts[idx2] = np.sum(spike_clusters == cluster_id2)
  284. if num_spikes_in_cluster > max_spikes_for_cluster:
  285. relative_counts = spike_counts / num_spikes_in_cluster * max_spikes_for_cluster
  286. else:
  287. relative_counts = spike_counts
  289. all_pcs = np.zeros((0, pc_features.shape[1], channels_to_use.size))
  290. all_labels = np.zeros((0,))
  292. for idx2, cluster_id2 in enumerate(units_in_range):
  294. subsample = int(relative_counts[idx2])
  296. pcs = get_unit_pcs(cluster_id2, spike_clusters, spike_templates,
  297. pc_feature_ind, pc_features, channels_to_use,
  298. subsample)
  300. if pcs is not None and len(pcs.shape) == 3:
  302. labels = np.ones((pcs.shape[0],)) * cluster_id2
  304. all_pcs = np.concatenate((all_pcs, pcs),0)
  305. all_labels = np.concatenate((all_labels, labels),0)
  307. all_pcs = np.reshape(all_pcs, (all_pcs.shape[0], pc_features.shape[1]*channels_to_use.size))
  308. if ((all_pcs.shape[0] > 10)
  309. and not (all_labels == cluster_id).all() # Not all labels are this cluster
  310. and (sum(all_labels == cluster_id) > 20) # No fewer than 20 spikes in this cluster
  311. and (len(channels_to_use) > 0)):
  312. isolation_distance, l_ratio = mahalanobis_metrics(all_pcs, all_labels, cluster_id)
  314. d_prime = lda_metrics(all_pcs, all_labels, cluster_id)
  316. nn_hit_rate, nn_miss_rate = nearest_neighbors_metrics(all_pcs, all_labels,
  317. cluster_id,
  318. max_spikes_for_nn,
  319. n_neighbors)
  320. else: # Too few spikes or cluster doesnt exist
  321. isolation_distance = np.nan
  322. d_prime = np.nan
  323. nn_miss_rate = np.nan
  324. nn_hit_rate = np.nan
  325. l_ratio = np.nan
  326. return isolation_distance, d_prime, nn_miss_rate, nn_hit_rate, l_ratio
  329. def calculate_pc_metrics(spike_clusters,
  330. spike_templates,
  331. total_units,
  332. pc_features,
  333. pc_feature_ind,
  334. num_channels_to_compare,
  335. max_spikes_for_cluster,
  336. max_spikes_for_nn,
  337. n_neighbors,
  338. do_parallel=True):
  339. """
  341. :param spike_clusters:
  342. :param total_units:
  343. :param pc_features:
  344. :param pc_feature_ind:
  345. :param num_channels_to_compare:
  346. :param max_spikes_for_cluster:
  347. :param max_spikes_for_nn:
  348. :param n_neighbors:
  349. :return:
  350. """
  352. assert (num_channels_to_compare % 2 == 1)
  353. half_spread = int((num_channels_to_compare - 1) / 2)
  355. cluster_ids = np.unique(spike_clusters)
  356. template_ids = np.unique(spike_templates)
  358. template_peak_channels = np.zeros((len(template_ids),), dtype='uint16')
  359. cluster_peak_channels = np.zeros((len(cluster_ids),), dtype='uint16')
  361. for idx, template_id in enumerate(template_ids):
  362. for_template = np.squeeze(spike_templates == template_id)
  363. pc_max = np.argmax(np.mean(pc_features[for_template, 0, :], 0))
  364. template_peak_channels[idx] = pc_feature_ind[template_id, pc_max]
  366. for idx, cluster_id in enumerate(cluster_ids):
  367. for_unit = np.squeeze(spike_clusters == cluster_id)
  368. templates_for_unit = np.unique(spike_templates[for_unit])
  369. template_positions = np.where(np.isin(template_ids, templates_for_unit))[0]
  370. cluster_peak_channels[idx] = np.median(template_peak_channels[template_positions])
  372. # Loop over clusters:
  373. if do_parallel:
  374. from joblib import Parallel, delayed
  375. meas = Parallel(n_jobs=-1, verbose=3)( # -1 means use all cores
  376. delayed(calculate_pc_metrics_one_cluster) # Function
  377. (cluster_peak_channels, idx, cluster_id, cluster_ids,
  378. half_spread, pc_features, pc_feature_ind,
  379. spike_clusters, spike_templates,
  380. max_spikes_for_cluster, max_spikes_for_nn, n_neighbors
  381. )
  382. for idx, cluster_id in enumerate(cluster_ids)) # Loop
  383. else:
  384. from tqdm import tqdm
  385. meas = []
  386. for idx, cluster_id in tqdm(enumerate(cluster_ids), total=cluster_ids.max(), desc='PC metrics'): # Loop
  387. meas.append(calculate_pc_metrics_one_cluster( # Function
  388. cluster_peak_channels, idx, cluster_id, cluster_ids,
  389. half_spread, pc_features, pc_feature_ind,
  390. spike_clusters, spike_templates,
  391. max_spikes_for_cluster, max_spikes_for_nn, n_neighbors))
  393. # Unpack:
  394. isolation_distances = []
  395. l_ratios = []
  396. d_primes = []
  397. nn_hit_rates = []
  398. nn_miss_rates = []
  399. for mea in meas:
  400. isolation_distance, d_prime, nn_miss_rate, nn_hit_rate, l_ratio = mea
  401. isolation_distances.append(isolation_distance)
  402. d_primes.append(d_prime)
  403. nn_miss_rates.append(nn_miss_rate)
  404. nn_hit_rates.append(nn_hit_rate)
  405. l_ratios.append(l_ratio)
  407. return (np.array(isolation_distances), np.array(l_ratios), np.array(d_primes),
  408. np.array(nn_hit_rates), np.array(nn_miss_rates))
  411. def calculate_silhouette_score(spike_clusters,
  412. spike_templates,
  413. total_units,
  414. pc_features,
  415. pc_feature_ind,
  416. total_spikes,
  417. do_parallel=True):
  419. def score_inner_loop(i, cluster_ids):
  420. """
  421. Helper to loop over cluster_ids in one dimension. We dont want to loop over both dimensions in parallel-
  422. that will create too much worker overhead
  423. Args:
  424. i: index of first dimension
  425. cluster_ids: iterable of cluster ids
  427. Returns: scores for dimension j
  429. """
  430. scores_1d = []
  431. for j in cluster_ids:
  432. if j > i:
  433. inds = np.in1d(cluster_labels, np.array([i, j]))
  434. X = all_pcs[inds, :]
  435. labels = cluster_labels[inds]
  437. # len(np.unique(labels))=1 Can happen if total_spikes is low:
  438. if (len(labels) > 2) and (len(np.unique(labels)) > 1):
  439. scores_1d.append(silhouette_score(X, labels))
  440. else:
  441. scores_1d.append(np.nan)
  442. else:
  443. scores_1d.append(np.nan)
  444. return scores_1d
  446. cluster_ids = np.unique(spike_clusters)
  448. random_spike_inds = np.random.permutation(spike_clusters.size)
  449. random_spike_inds = random_spike_inds[:total_spikes]
  450. num_pc_features = pc_features.shape[1]
  451. num_channels = np.max(pc_feature_ind) + 1
  453. all_pcs = np.zeros((total_spikes, num_channels * num_pc_features))
  455. for idx, i in enumerate(random_spike_inds):
  457. unit_id = spike_templates[i]
  458. channels = pc_feature_ind[unit_id,:]
  460. for j in range(0,num_pc_features):
  461. all_pcs[idx, channels + num_channels * j] = pc_features[i,j,:]
  463. cluster_labels = np.squeeze(spike_clusters[random_spike_inds])
  465. SS = np.empty((total_units, total_units))
  466. SS[:] = np.nan
  469. # Build lists
  470. if do_parallel:
  471. from joblib import Parallel, delayed
  472. scores = Parallel(n_jobs=-1, verbose=2)(delayed(score_inner_loop)(i, cluster_ids) for i in cluster_ids)
  473. else:
  474. scores = [score_inner_loop(i, cluster_ids) for i in cluster_ids]
  476. # Fill the 2d array
  477. for i, col_score in enumerate(scores):
  478. for j, one_score in enumerate(col_score):
  479. SS[i, j] = one_score
  481. with warnings.catch_warnings():
  482. warnings.simplefilter("ignore")
  483. a = np.nanmin(SS, 0)
  484. b = np.nanmin(SS, 1)
  486. return np.array([np.nanmin([a,b]) for a, b in zip(a,b)])
  489. def calculate_drift_metrics(spike_times,
  490. spike_clusters,
  491. spike_templates,
  492. total_units,
  493. pc_features,
  494. pc_feature_ind,
  495. interval_length,
  496. min_spikes_per_interval,
  497. do_parallel=True):
  498. def calc_one_cluster(cluster_id):
  499. """
  500. Helper to calculate drift for one cluster
  501. Args:
  502. cluster_id:
  504. Returns:
  505. max_drift, cumulative_drift
  506. """
  507. in_cluster = spike_clusters == cluster_id
  508. times_for_cluster = spike_times[in_cluster]
  509. depths_for_cluster = depths[in_cluster]
  511. median_depths = []
  513. for t1, t2 in zip(interval_starts, interval_ends):
  515. in_range = (times_for_cluster > t1) * (times_for_cluster < t2)
  517. if np.sum(in_range) >= min_spikes_per_interval:
  518. median_depths.append(np.median(depths_for_cluster[in_range]))
  519. else:
  520. median_depths.append(np.nan)
  522. median_depths = np.array(median_depths)
  523. max_drift = np.around(np.nanmax(median_depths) - np.nanmin(median_depths), 2)
  524. cumulative_drift = np.around(np.nansum(np.abs(np.diff(median_depths))), 2)
  525. return max_drift, cumulative_drift
  528. max_drifts = []
  529. cumulative_drifts = []
  531. depths = get_spike_depths(spike_templates, pc_features, pc_feature_ind)
  533. interval_starts = np.arange(np.min(spike_times), np.max(spike_times), interval_length)
  534. interval_ends = interval_starts + interval_length
  536. cluster_ids = np.unique(spike_clusters)
  538. if do_parallel:
  539. from joblib import Parallel, delayed
  540. meas = Parallel(n_jobs=-1, verbose=2)(delayed(calc_one_cluster)(cluster_id)
  541. for cluster_id in cluster_ids)
  542. else:
  543. meas = [calc_one_cluster(cluster_id) for cluster_id in cluster_ids]
  545. for max_drift, cumulative_drift in meas:
  546. max_drifts.append(max_drift)
  547. cumulative_drifts.append(cumulative_drift)
  548. return np.array(max_drifts), np.array(cumulative_drifts)
  552. # ==========================================================
  554. # IMPLEMENTATION OF ACTUAL METRICS:
  556. # ==========================================================
  559. def isi_violations(spike_train, min_time, max_time, isi_threshold, min_isi=0):
  560. """Calculate ISI violations for a spike train.
  562. Based on metric described in Hill et al. (2011) J Neurosci 31: 8699-8705
  564. modified by Dan Denman from cortex-lab/sortingQuality GitHub by Nick Steinmetz
  566. Inputs:
  567. -------
  568. spike_train : array of spike times
  569. min_time : minimum time for potential spikes
  570. max_time : maximum time for potential spikes
  571. isi_threshold : threshold for isi violation
  572. min_isi : threshold for duplicate spikes
  574. Outputs:
  575. --------
  576. fpRate : rate of contaminating spikes as a fraction of overall rate
  577. A perfect unit has a fpRate = 0
  578. A unit with some contamination has a fpRate < 0.5
  579. A unit with lots of contamination has a fpRate > 1.0
  580. num_violations : total number of violations
  582. """
  584. duplicate_spikes = np.where(np.diff(spike_train) <= min_isi)[0]
  586. spike_train = np.delete(spike_train, duplicate_spikes + 1)
  587. isis = np.diff(spike_train)
  589. num_spikes = len(spike_train)
  590. num_violations = sum(isis < isi_threshold)
  591. violation_time = 2*num_spikes*(isi_threshold - min_isi)
  592. total_rate = firing_rate(spike_train, min_time, max_time)
  593. violation_rate = num_violations/violation_time
  594. fpRate = violation_rate/total_rate
  596. return fpRate, num_violations
  599. def isi_violations_corrected(spike_train, min_time, max_time, isi_threshold, min_isi=0):
  600. """Calculate ISI violations for a spike train (with bias correction).
  602. This function was updated in September 2023 by Nick Steinmetz to correct
  603. two problems with the original implementation (text copied from
  604. https://github.com/cortex-lab/sortingQuality repo):
  606. 1) The approximation previously used, which was chosen to avoid getting
  607. imaginary results, wasn't accurate to the Hill et al paper on which this
  608. method was based, nor was it accurate to the correct solution to the problem.
  610. 2) Hill et al also did not have the correct solution to the problem. The Hill
  611. paper used an expression derived from an earlier work (Meunier et al
  612. 2003) which had assumed a special case: the "contamination" was itself
  613. only generated by a single neuron and therefore the contaminating spikes
  614. themselves had a refractory period. If instead the contaminating spikes
  615. are generated from a real Poisson process (as in the case of electrical
  616. noise or many nearby neurons generating the contamination), then the
  617. correct expression is different, as now calculated here. This expression
  618. is given in Llobet et al. bioRxiv 2022:
  620. https://www.biorxiv.org/content/10.1101/2022.02.08.479192v1.full.pdf
  622. In practice, the three methods (the real Hill equation, the original
  623. isi_violations calculation, and the correct equation implemented below)
  624. return almost identical values for contamination less than ~20%. They
  625. diverge strongly for 30% or more. For contamination levels above 50%
  626. (based on the original calculation), the corrected version is undefined
  627. due to a negative square root.
  629. Inputs:
  630. -------
  631. spike_train : array of spike times
  632. min_time : minimum time for potential spikes
  633. max_time : maximum time for potential spikes
  634. isi_threshold : threshold for isi violation
  635. min_isi : threshold for duplicate spikes
  637. Outputs:
  638. --------
  639. fpRate : rate of contaminating spikes as a fraction of overall rate
  640. A perfect unit has a fpRate = 0
  641. A unit with some contamination has a fpRate < 0.5
  642. A unit with lots of contamination has a fpRate > 1.0
  643. num_violations : total number of violations
  645. """
  647. duplicate_spikes = np.where(np.diff(spike_train) <= min_isi)[0]
  649. spike_train = np.delete(spike_train, duplicate_spikes + 1)
  650. isis = np.diff(spike_train)
  651. duration = max_time - min_time
  653. num_spikes = len(spike_train)
  654. num_violations = sum(isis < isi_threshold)
  655. fpRate = 1 - np.sqrt(1 - num_violations * duration /
  656. (pow(num_spikes, 2) * (isi_threshold - min_isi)))
  658. return fpRate, num_violations
  661. def presence_ratio(spike_train, min_time, max_time, num_bins=100):
  662. """Calculate fraction of time the unit is present within an epoch.
  664. Inputs:
  665. -------
  666. spike_train : array of spike times
  667. min_time : minimum time for potential spikes
  668. max_time : maximum time for potential spikes
  670. Outputs:
  671. --------
  672. presence_ratio : fraction of time bins in which this unit is spiking
  674. """
  676. h, b = np.histogram(spike_train, np.linspace(min_time, max_time, num_bins))
  678. return np.sum(h > 0) / (num_bins - 1)
  681. def firing_rate(spike_train, min_time = None, max_time = None):
  682. """Calculate firing rate for a spike train.
  684. If no temporal bounds are specified, the first and last spike time are used.
  686. Inputs:
  687. -------
  688. spike_train : numpy.ndarray
  689. Array of spike times in seconds
  690. min_time : float
  691. Time of first possible spike (optional)
  692. max_time : float
  693. Time of last possible spike (optional)
  695. Outputs:
  696. --------
  697. fr : float
  698. Firing rate in Hz
  700. """
  702. if min_time is not None and max_time is not None:
  703. duration = max_time - min_time
  704. else:
  705. duration = np.max(spike_train) - np.min(spike_train)
  707. fr = spike_train.size / duration
  709. return fr
  712. def amplitude_cutoff(amplitudes, num_histogram_bins = 500, histogram_smoothing_value = 3):
  714. """ Calculate approximate fraction of spikes missing from a distribution of amplitudes
  716. Assumes the amplitude histogram is symmetric (not valid in the presence of drift)
  718. Inspired by metric described in Hill et al. (2011) J Neurosci 31: 8699-8705
  720. Input:
  721. ------
  722. amplitudes : numpy.ndarray
  723. Array of amplitudes (don't need to be in physical units)
  725. Output:
  726. -------
  727. fraction_missing : float
  728. Fraction of missing spikes (0-0.5)
  729. If more than 50% of spikes are missing, an accurate estimate isn't possible
  731. """
  734. h,b = np.histogram(amplitudes, num_histogram_bins, density=True)
  736. pdf = gaussian_filter1d(h,histogram_smoothing_value)
  737. support = b[:-1]
  739. peak_index = np.argmax(pdf)
  740. G = np.argmin(np.abs(pdf[peak_index:] - pdf[0])) + peak_index
  742. bin_size = np.mean(np.diff(support))
  743. fraction_missing = np.sum(pdf[G:])*bin_size
  745. fraction_missing = np.min([fraction_missing, 0.5])
  747. return fraction_missing
  750. def mahalanobis_metrics(all_pcs, all_labels, this_unit_id):
  752. """ Calculates isolation distance and L-ratio (metrics computed from Mahalanobis distance)
  754. Based on metrics described in Schmitzer-Torbert et al. (2005) Neurosci 131: 1-11
  756. Inputs:
  757. -------
  758. all_pcs : numpy.ndarray (num_spikes x PCs)
  759. 2D array of PCs for all spikes
  760. all_labels : numpy.ndarray (num_spikes x 0)
  761. 1D array of cluster labels for all spikes
  762. this_unit_id : Int
  763. number corresponding to unit for which these metrics will be calculated
  765. Outputs:
  766. --------
  767. isolation_distance : float
  768. Isolation distance of this unit
  769. l_ratio : float
  770. L-ratio for this unit
  772. """
  774. pcs_for_this_unit = all_pcs[all_labels == this_unit_id,:]
  775. pcs_for_other_units = all_pcs[all_labels != this_unit_id, :]
  777. mean_value = np.expand_dims(np.mean(pcs_for_this_unit,0),0)
  779. try:
  780. VI = np.linalg.inv(np.cov(pcs_for_this_unit.T))
  781. except np.linalg.linalg.LinAlgError: # case of singular matrix
  782. return np.nan, np.nan
  784. mahalanobis_other = np.sort(cdist(mean_value,
  785. pcs_for_other_units,
  786. 'mahalanobis', VI = VI)[0])
  788. mahalanobis_self = np.sort(cdist(mean_value,
  789. pcs_for_this_unit,
  790. 'mahalanobis', VI = VI)[0])
  792. n = np.min([pcs_for_this_unit.shape[0], pcs_for_other_units.shape[0]]) # number of spikes
  794. if n >= 2:
  796. dof = pcs_for_this_unit.shape[1] # number of features
  798. l_ratio = np.sum(1 - chi2.cdf(pow(mahalanobis_other,2), dof)) / \
  799. mahalanobis_self.shape[0] # normalize by size of cluster, not number of other spikes
  800. isolation_distance = pow(mahalanobis_other[n-1],2)
  802. else:
  803. l_ratio = np.nan
  804. isolation_distance = np.nan
  806. return isolation_distance, l_ratio
  811. def lda_metrics(all_pcs, all_labels, this_unit_id):
  813. """ Calculates d-prime based on Linear Discriminant Analysis
  815. Based on metric described in Hill et al. (2011) J Neurosci 31: 8699-8705
  817. Inputs:
  818. -------
  819. all_pcs : numpy.ndarray (num_spikes x PCs)
  820. 2D array of PCs for all spikes
  821. all_labels : numpy.ndarray (num_spikes x 0)
  822. 1D array of cluster labels for all spikes
  823. this_unit_id : Int
  824. number corresponding to unit for which these metrics will be calculated
  826. Outputs:
  827. --------
  828. d_prime : float
  829. Isolation distance of this unit
  830. l_ratio : float
  831. L-ratio for this unit
  833. """
  835. X = all_pcs
  837. y = np.zeros((X.shape[0],),dtype='bool')
  838. y[all_labels == this_unit_id] = True
  840. lda = LDA(n_components=1)
  842. X_flda = lda.fit_transform(X, y)
  844. flda_this_cluster = X_flda[np.where(y)[0]]
  845. flda_other_cluster = X_flda[np.where(np.invert(y))[0]]
  847. d_prime = (np.mean(flda_this_cluster) - np.mean(flda_other_cluster))/np.sqrt(0.5*(np.std(flda_this_cluster)**2+np.std(flda_other_cluster)**2))
  849. return d_prime
  853. def nearest_neighbors_metrics(all_pcs, all_labels, this_unit_id, max_spikes_for_nn, n_neighbors):
  855. """ Calculates unit contamination based on NearestNeighbors search in PCA space
  857. Based on metrics described in Chung, Magland et al. (2017) Neuron 95: 1381-1394
  859. Inputs:
  860. -------
  861. all_pcs : numpy.ndarray (num_spikes x PCs)
  862. 2D array of PCs for all spikes
  863. all_labels : numpy.ndarray (num_spikes x 0)
  864. 1D array of cluster labels for all spikes
  865. this_unit_id : Int
  866. number corresponding to unit for which these metrics will be calculated
  867. max_spikes_for_nn : Int
  868. number of spikes to use (calculation can be very slow when this number is >20000)
  869. n_neighbors : Int
  870. number of neighbors to use
  872. Outputs:
  873. --------
  874. hit_rate : float
  875. Fraction of neighbors for target cluster that are also in target cluster
  876. miss_rate : float
  877. Fraction of neighbors outside target cluster that are in target cluster
  879. """
  881. total_spikes = all_pcs.shape[0]
  882. ratio = max_spikes_for_nn / total_spikes
  883. this_unit = all_labels == this_unit_id
  885. X = np.concatenate((all_pcs[this_unit,:], all_pcs[np.invert(this_unit),:]),0)
  887. n = np.sum(this_unit)
  889. if ratio < 1:
  890. inds = np.arange(0,X.shape[0]-1,1/ratio).astype('int')
  891. X = X[inds,:]
  892. n = int(n * ratio)
  895. nbrs = NearestNeighbors(n_neighbors=n_neighbors, algorithm='ball_tree').fit(X)
  896. distances, indices = nbrs.kneighbors(X)
  898. this_cluster_inds = np.arange(n)
  900. this_cluster_nearest = indices[:n,1:].flatten()
  901. other_cluster_nearest = indices[n:,1:].flatten()
  903. hit_rate = np.mean(this_cluster_nearest < n)
  904. miss_rate = np.mean(other_cluster_nearest < n)
  906. return hit_rate, miss_rate
  908. # ==========================================================
  910. # HELPER FUNCTIONS:
  912. # ==========================================================
  914. def features_intersect(pc_feature_ind, these_channels):
  915. """
  916. # Take only the channels that have calculated features out of the ones we are interested in:
  917. # This should reduce the occurence of 'except IndexError' below
  919. Args:
  920. pc_feature_ind
  921. these_channels: channels_to_use or units_for_channel
  923. Returns:
  924. channels_to_use: intersect of what's available in PCs and what's needed
  925. """
  926. intersect = set(pc_feature_ind[these_channels[0], :]) # Initialize
  927. for cluster_id2 in these_channels:
  928. # Make a running intersect of what is available and what is needed
  929. intersect = intersect & set(pc_feature_ind[cluster_id2, :])
  930. return np.array(list(intersect))
  933. def get_unit_pcs(unit_id,
  934. spike_clusters,
  935. spike_templates,
  936. pc_feature_ind,
  937. pc_features,
  938. channels_to_use,
  939. subsample):
  941. """ Return PC features for one unit
  943. Inputs:
  944. -------
  945. unit_id : Int
  946. ID for this unit
  947. spike_clusters : np.ndarray
  948. Cluster labels for each spike
  949. spike_templates : np.ndarry
  950. Template labels for each spike
  951. pc_feature_ind : np.ndarray
  952. Channels used for PC calculation for each unit
  953. pc_features : np.ndarray
  954. Array of all PC features
  955. channels_to_use : np.ndarray
  956. Channels to use for calculating metrics
  957. subsample : Int
  958. maximum number of spikes to return
  960. Output:
  961. -------
  962. unit_PCs : numpy.ndarray (float)
  963. PCs for one unit (num_spikes x num_PCs x num_channels)
  965. """
  968. inds_for_unit = np.where(spike_clusters == unit_id)[0]
  970. spikes_to_use = np.random.permutation(inds_for_unit)[:subsample]
  972. unique_template_ids = np.unique(spike_templates[spikes_to_use])
  974. unit_PCs = []
  976. for template_id in unique_template_ids:
  978. index_mask = spikes_to_use[np.squeeze(spike_templates[spikes_to_use]) == template_id]
  979. these_inds = pc_feature_ind[template_id, :]
  981. pc_array = []
  983. for i in channels_to_use:
  985. if np.isin(i, these_inds):
  986. channel_index = np.argwhere(these_inds == i)[0][0]
  987. pc_array.append(pc_features[index_mask, :, channel_index])
  988. else:
  989. return None
  991. unit_PCs.append(np.stack(pc_array, axis=-1))
  993. if len(unit_PCs) > 0:
  995. return np.concatenate(unit_PCs)
  996. else:
  997. return None

metrics.py at commit 9199927, under BSD-2-Clause · at the source

Overview

Authors: Julian Keil1,2,3, Victor Hernandez-Urbina2, Chrystalleni Vassiliou4,5, Camin Dean4,5,6,7, Dietmar Schmitz4,5,6,7,8,9, Jens Kremkow7,10, Jérémie Sibille4,5,6,7,8,9
ORCID iDs: Dietmar Schmitz
  1. Department of Cognitive Science, University of Potsdam, 14469 Potsdam, Germany
  2. Nuuron GmbH Berlin, Berlin, Germany
  3. Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität Berlin, Clinical Neurotechnology Laboratory, Department of Psychiatry and Neurosciences, 10117 Berlin, Germany
  4. German Center for Neurodegenerative Diseases (DZNE) Berlin, 10117 Berlin, Germany
  5. Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität Berlin, Neuroscience Research Center, 10117 Berlin, Germany
  6. Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität Berlin, Einstein Center for Neurosciences, 10117 Berlin, Germany
  7. Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität Berlin, Bernstein Center for Computational Neuroscience, 10115 Berlin, Germany
  8. Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität Berlin, NeuroCure Cluster of Excellence, 10117 Berlin, Germany
  9. Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität zu Berlin, Institute of Cell and Neurobiology, 10115 Berlin, Germany
  10. Institute of Biology, Otto von Guericke University Magdeburg, 39120 Magdeburg, Germany
Institutions: University of Potsdam (Germany); Humboldt-Universität zu Berlin (Germany); Charité - Universitätsmedizin Berlin (Germany); German Center for Neurodegenerative Diseases (Germany); Bernstein Center for Computational Neuroscience Berlin (Germany); NeuroCure Cluster of Excellence (Germany); Otto-von-Guericke-Universität Magdeburg (Germany)
Journal: Scientific reports, volume 16, issue 1, article 15228
Dates: received 26 January 2026; accepted 4 May 2026; published online 17 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41598-026-52253-9 · PMID 42144416 · PMCID PMC13181131 · OpenAlex W4406545532
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.01.14.632930
Status: code verified
Categories: mouse (organism), systems (subfield)
Methods: Spectral & time-frequency, Statistics, Preprocessing, Evoked potentials, Single-unit activity, calcium imaging
Keywords: Biological techniques, Neuroscience
MeSH: Photic Stimulation*, Visual Cortex*, Animals, Mice (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Charité - Universitätsmedizin Berlin
Citations: cited by 1 paper (Europe PMC); 43 references in the paper

Abstract

Today, flickering full-field visual stimulation is used to increase neuronal oscillations for a variety of research or therapeutic purposes. We propose spatially organized sequential visual flickering stimulation as a newer tool to circumvent the intrinsic low-pass filter nature of the vertebrate visual system, in order to increase the power of high frequency oscillations in the visual system. We show that spatially organized visual flickering can increase power in high frequencies (100 to 190 Hz) in the visual cortex of mice. Consequently, spatially organized sequential sensory stimulation should be regarded as a putative new way leading to power increases in high frequency domains.

Supplementary Information: The online version contains supplementary material available at 10.1038/s41598-026-52253-9.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

AllenInstitute/ecephys_spike_sorting

License: BSD-2-Clause
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 919992748a5324724ba87169ecdcf9eb6e3b9973, 18 August 2024
Languages: Python (77), C++ (19), C/C++ (10), JavaScript (5), Shell (2), MATLAB (1)
Size: 216 files, 114 scripts
Software Heritage: not archived
Found in: the text, “Hardware, software”
Holds: README, license file, environment (Pipfile, Pipfile.lock, setup.cfg, setup.py), tests, continuous integration, documentation
Not found: CITATION.cff
Tools: NumPy (34 files), pandas (11 files), SciPy (9 files), Matplotlib (8 files), scikit-learn (4 files), xarray (2 files), h5py (1 file), Kilosort (1 file), Pillow (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
116 files
At the source: github.com/AllenInstitute/ecephys_spike_sorting

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 114 scripts, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Data and scripts will be provided freely upon reasonable request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 2 keywords, 4 MeSH terms, 1 funder, 36 references.

Cite

This paper

Keil, J., Hernandez-Urbina, V., Vassiliou, C., Dean, C., Schmitz, D., Kremkow, J., & Sibille, J. (2026). Sequential visual stimuli increase high frequency power in the visual cortex. Scientific reports, 16(1), 15228. https://doi.org/10.1038/s41598-026-52253-9

BibTeX

@article{keil2026sequential,
author = {Keil, Julian and Hernandez-Urbina, Victor and Vassiliou, Chrystalleni and Dean, Camin and Schmitz, Dietmar and Kremkow, Jens and Sibille, Jérémie},
title = {{Sequential visual stimuli increase high frequency power in the visual cortex}},
journal = {Scientific reports},
year = {2026},
month = may,
volume = {16},
number = {1},
pages = {15228},
publisher = {Nature Publishing Group},
issn = {2045-2322},
doi = {10.1038/s41598-026-52253-9},
url = {https://doi.org/10.1038/s41598-026-52253-9},
pmid = {42144416},
pmcid = {PMC13181131}
}

RIS

TY - JOUR
AU - Keil, Julian
AU - Hernandez-Urbina, Victor
AU - Vassiliou, Chrystalleni
AU - Dean, Camin
AU - Schmitz, Dietmar
AU - Kremkow, Jens
AU - Sibille, Jérémie
TI - Sequential visual stimuli increase high frequency power in the visual cortex
T2 - Scientific reports
J2 - Sci Rep
PY - 2026
DA - 2026/05/17
VL - 16
IS - 1
SP - 15228
SN - 2045-2322
PB - Nature Publishing Group
DO - 10.1038/s41598-026-52253-9
UR - https://doi.org/10.1038/s41598-026-52253-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41598-026-52253-9",
"type": "article-journal",
"title": "Sequential visual stimuli increase high frequency power in the visual cortex",
"container-title": "Scientific reports",
"author": [
{
"family": "Keil",
"given": "Julian"
},
{
"family": "Hernandez-Urbina",
"given": "Victor"
},
{
"family": "Vassiliou",
"given": "Chrystalleni"
},
{
"family": "Dean",
"given": "Camin"
},
{
"family": "Schmitz",
"given": "Dietmar"
},
{
"family": "Kremkow",
"given": "Jens"
},
{
"family": "Sibille",
"given": "Jérémie"
}
],
"container-title-short": "Sci Rep",
"volume": "16",
"issue": "1",
"page": "15228",
"DOI": "10.1038/s41598-026-52253-9",
"PMID": "42144416",
"PMCID": "PMC13181131",
"ISSN": "2045-2322",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41598-026-52253-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
17
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: Kilosort, xarray, h5py, 6 other tools, systems, mouse, 1 reference
[2] doi:10.1162/imag.a.1229 [code]
40 Hz audiovisual stimulation improves sustained attention and related brain oscillations.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: pandas, SciPy, Matplotlib, 1 other tool, 7 references
[3] doi:10.1038/s41467-026-76581-6 [code]
Thalamocortical bursts encode reward contingencies and drive associative learning.
Journal: Nature communications
In common: Kilosort, xarray, h5py, 6 other tools, systems, mouse
[4] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: Kilosort, h5py, Pillow, 5 other tools, mouse, 2 references
[5] doi:10.1038/s41593-026-02258-4 [code]
Laminar organization of cellular microcircuits modulating human interictal epileptiform discharges.
Journal: Nature neuroscience
In common: Kilosort, scikit-learn, pandas, 3 other tools, systems, 3 references
[6] doi:10.1016/j.patter.2026.101590 [code]
Density-based longitudinal neuron tracking in high-density electrophysiological recordings.
Journal: Patterns (New York, N.Y.)
In common: Kilosort, h5py, Pillow, 5 other tools, 1 reference
[7] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: Kilosort, h5py, Pillow, 5 other tools, systems, mouse
[8] doi:10.1038/s41467-026-72057-9 [code]
Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: xarray, h5py, Pillow, 5 other tools, systems
[9] doi:10.1038/s41586-026-10331-y [code]
Active dissociation of intracortical spiking and high gamma activity.
Journal: Nature
In common: h5py, Pillow, scikit-learn, 4 other tools, 2 references
[10] doi:10.1098/rstb.2024.0461 [code]
Shallow recurrent decoders for neural and behavioural dynamics.
Journal: Philosophical transactions of the Royal Society of London. Series B, Biological sciences
In common: xarray, h5py, Pillow, 5 other tools, mouse

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.