OSCR

A flexible quality metric for electrophysiological recordings across brain regions and species

Code ↔ Paper

16 matches 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 16 matches · 5 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › False acceptance rate in the Sliding RP metric ↔ matlab/computeMatrix.m, lines 2–59 · score 0.79 · Markov chain, error rate, expected violations, passage, family, Sliding RP
  2. [2] § Methods › False acceptance rate in the Sliding RP metric ↔ python/slidingRP/metrics.py, lines 468–522 · score 0.76 · Markov chain, passage probability, error rate, family, Sliding RP, Poisson
  3. [3] § Methods › Estimating RP duration for comparison across regions and species (Fig. 1) ↔ roth-et-al-2026/rpQuantification/estimate_refractory_period.m, lines 45–55 · score 0.70 · peak detection, 0–0.5, narrow, outside, refractory period, filtered
  4. [4] § Methods › Algorithm Description › Step 4: Computing Confidence Scores ↔ matlab/computeViol.m, the whole file · a weak match · score 0.69 · fewer violations, violations Ve, RP violations, cumulative, Scores, probability
  5. [5] § Methods › Estimating RP duration for comparison across regions and species (Fig. 1) ↔ python/slidingRP/metrics.py, lines 34–126 · score 0.66 · 0–0.5, filtered ACG, global, outside, peak, milliseconds
  6. [6] § Methods › Algorithm Description › Step 2: Computing Expected Violations ↔ matlab/computeViol.m, the whole file · a weak match · score 0.65 · neuron spikes, contamination proportion, RP violations, spike train, hypothesis, contaminating spikes
  7. [7] § Methods › Estimating RP duration for comparison across regions and species (Fig. 1) ↔ roth-et-al-2026/rpQuantification/estimate_refractory_period.m, lines 57–79 · score 0.63 · valid peaks, refractory period, earliest, closest, filtered, ACG
  8. [8] § Results › A quality metric suitable for variable RP durations ↔ matlab/slidingRP_all.m, the whole file · a weak match · score 0.63 · maximum acceptable contamination, maximum confidence, Sliding RP metric, lowest, variable, matrix
  9. [9] § Methods › Estimating RP duration for comparison across regions and species (Fig. 1) ↔ roth-et-al-2026/rpQuantification/estimate_refractory_period.m, lines 92–136 · score 0.60 · fit parameters, midpoint, steepness, refractory period, ymax, ymin
  10. [10] § Results › Influence of Confidence parameter ↔ matlab/RPmetric_Classic.m, the whole file · a weak match · score 0.59 · statistical power, Hill Llobet, contamination threshold, matching, firing rate, Sliding RP
  11. [11] § Methods › Estimating RP duration for comparison across regions and species (Fig. 1) ↔ python/slidingRP/metrics.py, lines 34–126 · score 0.58 · filtered ACG, aborted, strictly, closest, baseline, peaks
  12. [12] § Results › A quality metric suitable for variable RP durations ↔ matlab/slidingRP.m, lines 2–68 · score 0.58 · maximum acceptable contamination, maximum confidence, Sliding RP, lowest, matrix, quality
  13. [13] § Methods › Algorithm Description › Step 1: Autocorrelogram Computation ↔ python/slidingRP/metrics.py, lines 361–398 · score 0.57 · spike pairs, nACG, histogram, lag, autocorrelogram, bin
  14. [14] § Results › Validation of the metric with simulated data ↔ matlab/slidingRP.m, lines 2–68 · score 0.56 · low firing rates, maximum acceptable, S3, spike trains, conservatively, Sliding RP
  15. [15] § Methods › Spike train simulations ↔ matlab/simulations/genST.m, the whole file · a weak match · score 0.55 · cumulative sum, exprnd, exponentially, spike train, vector, recording duration
  16. [16] § Methods › Estimating RP duration for comparison across regions and species (Fig. 1) ↔ roth-et-al-2026/rpQuantification/estimate_refractory_period.m, lines 36–43 · score 0.52 · median filter, smoothed, milliseconds, window, bin, ACG

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 · 776 lines · 32 KB · no license · 4 matches

  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Sun Jul 10 11:34:59 2022
  5. @author: Gaelle Chapuis ; Legacy code from Noam Roth commented below
  7. compute the metric for a single cluster (neuron) in a recording
  8. """
  9. import warnings
  10. from scipy.optimize import OptimizeWarning
  11. import numpy as np
  12. from scipy import stats
  13. from scipy.optimize import curve_fit
  14. import scipy
  17. def closest(lst, K):
  18. lst = np.asarray(lst)
  19. idx = (np.abs(lst - K)).argmin()
  20. return idx, lst[idx]
  23. def compute_timebins(acg, bin_size_secs):
  24. x = np.arange(len(acg))
  25. timeBins = x.dot(bin_size_secs)
  26. return timeBins
  29. def sigmoid(x, L, x0, k, b): # L is max value; x0 is the midpoint x value; k is the steepness, b is the baseline shift
  30. y = L / (1 + np.exp(-k * (x - x0))) + b
  31. return y
  34. def compute_rf(acg,
  35. min_sig=np.array([0.0, 0.0005]), # 0-0.5 ms
  36. t_medfilter=0.00083,
  37. bin_size_secs=1 / 30_000,
  38. timeBins=None,
  39. RPEstimateFromPercentageOfSlope = 0.05,
  40. fr_percentage=10 / 100):
  41. '''
  42. Compute the refractory period (RP) from an auto-correlogram (ACG) array, by :
  43. - filtering the ACG (median filter)
  44. - finding the ACG first local peak above a certain firing rate threshold
  45. - fitting a sigmoid to the first portion of the filtered ACG
  46. - The RP is defined as the time at which a certain percentage of the sigmoid fit is reached.
  47. :param acg: the autocorrelogram, 1D numpy array containing spike number in bins of size bin_size_secs
  48. :param min_sig: time window in second to compute the minimum of the sigmoid fit
  49. :param t_medfilter: window for the median filter, in seconds
  50. :param bin_size_secs: size of the acg bins, in seconds
  51. :param timeBins: time of the bins, re-computed if None
  52. :param RPEstimateFromPercentageOfSlope: portion of the sigmoid fit at which to take the time as RP; value from 0-1
  53. :param fr_percentage: percentage of the max-min firing rate, used to find the first peak.
  54. :return:
  55. estimatedRP: the estimated RP in milliseconds
  56. estimateIdx: the index of the RP (which corresponds to the bin index of the ACG)
  57. xSigmoid, ySigmoid: the sigmoid fit (x-values: times, y-values: curve values)
  59. '''
  60. # This is a hack: we do not want to return the fit on an ACG where the fit is poor
  61. # Treat a poor curve_fit as a failure, but only within this function
  62. # (catch_warnings restores the global filter state on exit).
  63. with warnings.catch_warnings():
  64. warnings.simplefilter("error", OptimizeWarning)
  66. if timeBins is None: # Compute timebins
  67. timeBins = compute_timebins(acg, bin_size_secs)
  69. # Median filter
  70. med_filt = scipy.ndimage.median_filter(acg, size=int(np.round(t_medfilter / bin_size_secs)))
  71. # Find all the peaks on this filtered trace
  72. peaks_idx = scipy.signal.find_peaks(med_filt)[0]
  74. # Compute value (max) in short time window near 0
  75. minSigmoid = np.max(med_filt[(timeBins > min_sig[0]) & (timeBins < min_sig[1])])
  76. # Note: could be using either the max() or mean()
  77. # max forces the peak to be outside of this time window
  79. # To find the peak, make sure it is strictly higher than this value
  80. # and above a baseline percentage firing rate
  81. fr_baseline = fr_percentage * (med_filt.max() - med_filt.min())
  82. peaks_possible = np.where((med_filt[peaks_idx] > minSigmoid) &
  83. (med_filt[peaks_idx] > fr_baseline))[0]
  85. # if no peak is found, abort and return NaN
  86. if len(peaks_possible) == 0:
  87. estimatedRP = np.nan
  88. estimateIdx = np.nan
  89. xSigmoid = np.array([])
  90. ySigmoid = np.array([])
  91. return estimatedRP, estimateIdx, xSigmoid, ySigmoid
  93. # Use first peak possible found (i.e. closest to 0 second on the ACG)
  94. peak_idx = peaks_idx[peaks_possible[0]]
  95. maxSigmoid = med_filt[peak_idx]
  97. # Truncate ACG and time bins according to max
  98. timeBins_fit = timeBins[0:peak_idx]
  99. acg_fit = med_filt[0:peak_idx]
  101. # fit the sigmoid with max and min fixed
  102. try:
  103. popt, pcov = curve_fit(lambda x, x0, k: sigmoid(x, maxSigmoid, x0, k, minSigmoid), timeBins_fit, acg_fit)
  104. fitParams = [maxSigmoid, popt[0], popt[1], minSigmoid]
  106. xSigmoid = timeBins
  107. ySigmoid = sigmoid(xSigmoid, *fitParams)
  109. # find RP
  110. estimateIdx, _ = closest(ySigmoid, RPEstimateFromPercentageOfSlope * (maxSigmoid - minSigmoid) + minSigmoid)
  112. # Compute the index of the first ACG bin with non-null firing rate
  113. first_index = np.where(acg != 0)[0][0]
  114. if estimateIdx < first_index: # If the estimate RP is BEFORE the first ACG bin
  115. estimateIdx = first_index # Replace
  117. estimatedRP = 1000 * xSigmoid[estimateIdx] # in ms
  118. if np.max(ySigmoid) - np.min(ySigmoid) < 1: # The fit is essentially flat
  119. raise OptimizeWarning
  120. except (OptimizeWarning, RuntimeWarning, RuntimeError): # This is in the case the bins are too few
  121. # print('fit error')
  122. estimatedRP = np.nan
  123. estimateIdx = np.nan
  124. xSigmoid = np.array([])
  125. ySigmoid = np.array([])
  126. return estimatedRP, estimateIdx, xSigmoid, ySigmoid
  129. def remove_lowrp_confmat(confMatrix, rp, rp_reject=0.0005):
  130. # We want to compute on the matrix only for RPs above a certain value
  131. # Remove those small RP values from the rp vector and conf matrix
  132. rp_idx_keep = rp > rp_reject
  133. rp = rp[rp_idx_keep]
  134. confMatrix = confMatrix[:, rp_idx_keep]
  135. return confMatrix, rp
  138. def confidence_contamin(confMatrix, cont, rp, cont_thresh=10.0, rp_reject = 0.0005):
  139. '''
  140. For a level of contamination contamin_level given (default 10%), find the smallest confidence
  141. value for which the minimum value of the contamination curve is equal or lower to the
  142. contamin_level.
  143. In practice, this is equivalent to finding the maximum value of the confidence of
  144. the confidence matrix at the contamination level row.
  146. Uses the output of the function computeMatrix()
  148. :param confMatrix: the confidence matrix (contamination x RP, values: confidence, ranging from 0-1)
  149. :param cont: contamination vector at which the confidence is computed (ranges by default from 0-35)
  150. :param rp: refractory period vector at which the confidence is computed
  151. :param cont_level: level of contamination searched for, default is 10% (0.1)
  152. :return:
  153. '''
  154. # Find index in cont vector where there is cont_thresh or closest (higher) value
  155. idx_cont = np.where(cont >= cont_thresh)[0][0]
  156. cont_thresh = cont[idx_cont] # Return actual level of contamination studied
  158. # We want to compute the curve of contamination only for RPs above a certain value
  159. # Remove those small RP values from the rp vector and conf matrix
  160. confMatrix, _ = remove_lowrp_confmat(confMatrix, rp, rp_reject=rp_reject)
  162. # At the contamination level studied, find the maximal value of confidence
  163. max_conf = np.max(confMatrix[idx_cont, :]) # Legacy name: 'maxConfidenceAt10Cont'
  164. return max_conf, idx_cont, cont_thresh
  167. def pass_slidingRP_confmat(confMatrix, cont, rp, conf_thresh=90, cont_thresh=10, rp_reject=0.0005):
  168. '''
  169. Given a confidence matrix, a confidence threshold (default 90%) and a contamination threshold (default=10),
  170. assess whether the unit passes the sliding RP metric
  172. Uses the output of the function computeMatrix()
  174. :param confMatrix:
  175. :param cont:
  176. :param rp:
  177. :param conf_thresh:
  178. :param cont_thresh:
  179. :param rp_reject:
  180. :return:
  181. '''
  182. # We want to compute the curve of contamination only for RPs above a certain value
  183. # Remove those small RP values from the rp vector and conf matrix
  184. confMatrix, rp = remove_lowrp_confmat(confMatrix, rp, rp_reject=rp_reject)
  186. # Find matrix indices that are above or equal to the confidence threshold
  187. a = np.where(confMatrix >= conf_thresh)
  188. if len(a[0]) > 0:
  189. # Find minimum contamination value for this conf threshold
  190. min_idx = np.min(a[0]) # Min on rows axis = contamination axis
  191. min_cont = cont[min_idx] # Legacy name: minContWith90Confidence
  192. # Find the smallest RP possible at the contamination level at the max confidence val
  193. minRP = np.argmax(confMatrix[min_idx, :])
  194. rp_min_val = rp[minRP] # Legacy name: timeOfLowestCont
  195. # rp_min_val is the tau_r at which the minimum contamination is confirmed,
  196. # i.e. an estimate of the unit's RP duration (manuscript Outputs, Step 6).
  198. # Check if this unit passes the sliding RP metric
  199. pass_cont_thresh = min_cont <= cont_thresh
  200. else:
  201. pass_cont_thresh = False
  202. min_cont = np.nan
  203. rp_min_val = np.nan
  205. return pass_cont_thresh, min_cont, rp_min_val # Legacy: value, minContWith90Confidence, timeOfLowestCont
  208. def slidingRP(spikeTimes, params=None, conf_thresh=90, cont_thresh=10, rp_reject=0.0005):
  209. """Compute the Sliding RP metric for a single cluster.
  211. Mirrors the MATLAB slidingRP.m. Pass options in the ``params`` dict
  212. (e.g. ``params={'recDur': 3600}``); recDur is the recording duration in
  213. seconds and defaults to ``max(spikeTimes)`` if not given (recommended to
  214. set explicitly). The ACG and expected-violation (Llobet) calculation are
  215. identical to the MATLAB implementation.
  217. The maximum confidence and pass/fail come from the per-tau_r confidence at
  218. the contamination threshold; the minimum confirmable contamination
  219. (``min_cont``) is computed analytically (closed-form, no contamination grid).
  220. For the full confidence matrix, call computeMatrix() directly.
  222. params options include: recDur, sampleRate, binSizeCorr, cont, correction
  223. (FWER multiple-comparisons; default off), and forcePass (IBL-specific; default
  224. off) — when True, units with zero <2 ms violations and firing_rate > 0.5 are
  225. flagged via pass_forced.
  227. Returns
  228. -------
  229. max_conf : float maximum confidence (%) that contamination < cont_thresh
  230. min_cont : float minimum contamination (%) confirmable at conf_thresh
  231. (continuous; NaN if > 35%)
  232. rp_min_val : float tau_r (s) at which min_cont is achieved
  233. n_spikes_below2 : int ACG count with ISI < 2 ms
  234. firing_rate : float n_spikes / recDur
  235. pass_cont_thresh : bool max_conf >= conf_thresh
  236. pass_forced : bool IBL force-pass flag (False unless params['forcePass'])
  238. Mapping to MATLAB slidingRP.m outputs (which differ in order):
  239. MATLAB [passTest, confidence, contamination, timeOfLowestCont, nViolShort, ...]
  240. Python pass_cont_thresh, max_conf, min_cont, rp_min_val, n_spikes_below2
  241. (MATLAB's confMatrix/cont/rp/nACG come from computeMatrix(); Python returns
  242. firing_rate and pass_forced instead. Values match; only order/names differ.)
  243. """
  244. params = dict(params) if params else {}
  245. sampleRate = params.setdefault('sampleRate', 30000)
  246. rpBinSize = params.setdefault('binSizeCorr', 1 / sampleRate)
  248. spikeTimes = np.asarray(spikeTimes, dtype=np.float64)
  250. if params.get('correction', False):
  251. # FWER-corrected variant (Fig. S3): derive the scalar metrics from the
  252. # corrected confidence matrix (grid-based; no analytical shortcut for the
  253. # corrected confidence). Expensive, non-default path.
  254. confMatrix, cont, rp, nACG, firing_rate = computeMatrix(spikeTimes, params)
  255. pass_cont_thresh, min_cont, rp_min_val = pass_slidingRP_confmat(
  256. confMatrix, cont, rp, conf_thresh, cont_thresh, rp_reject)
  257. max_conf, _, _ = confidence_contamin(confMatrix, cont, rp, cont_thresh, rp_reject)
  258. n_spikes_below2 = int(np.sum(nACG[0:np.where(rp > 0.002)[0][0] + 1]))
  259. pass_forced = params.get('forcePass', False) and (n_spikes_below2 == 0) \
  260. and (firing_rate > 0.5) and (not pass_cont_thresh)
  261. return max_conf, min_cont, rp_min_val, n_spikes_below2, firing_rate, \
  262. pass_cont_thresh, pass_forced
  264. n_spikes = spikeTimes.size
  265. recDur = params.get('recDur', None)
  266. if recDur is None:
  267. recDur = float(np.max(spikeTimes)) if n_spikes else 0.0
  269. rpEdges = np.arange(0, 10 / 1000, rpBinSize) # in s
  270. rp = rpEdges + rpBinSize / 2 # bin centres (s)
  271. refDur = rp + rpBinSize / 2 # right bin edge = tested tau_r
  272. nACG = computeACG(spikeTimes, rpBinSize, rp.size)
  273. obsViol = np.cumsum(nACG)
  274. firing_rate = n_spikes / recDur if recDur > 0 else 0.0
  276. testTimes = rp > rp_reject # exclude tau_r below rp_reject from pass/fail
  278. # Max confidence at the contamination threshold (identical to the matrix path)
  279. conf_at_thresh = 100 * computeViol(obsViol, firing_rate, n_spikes, refDur,
  280. cont_thresh / 100, recDur)
  281. max_conf = float(np.max(conf_at_thresh[testTimes])) if np.any(testTimes) else 0.0
  283. # Minimum confirmable contamination, computed analytically (continuous)
  284. min_cont, rp_min_val = compute_min_contamination(
  285. obsViol, n_spikes, refDur, rp, recDur, conf_thresh, rp_reject)
  287. pass_cont_thresh = bool(max_conf >= conf_thresh)
  289. n_spikes_below2 = int(np.sum(nACG[0:np.where(rp > 0.002)[0][0] + 1]))
  291. # IBL-specific: force-pass units with zero short-ISI violations and FR > 0.5
  292. # (tuned for IBL ~1 h recordings). Opt-in only (default off), matching the
  293. # correction path and the docstring.
  294. pass_forced = params.get('forcePass', False) and (n_spikes_below2 == 0) \
  295. and (firing_rate > 0.5) and (not pass_cont_thresh)
  297. return max_conf, min_cont, rp_min_val, \
  298. n_spikes_below2, firing_rate, \
  299. pass_cont_thresh, pass_forced
  302. def _slidingRP_worker(args):
  303. """Module-level worker so slidingRP_all can parallelise with a process pool
  304. (the callable and its args must be picklable)."""
  305. st, params, conf_thresh, cont_thresh, rp_reject = args
  306. return slidingRP(st, params=params, conf_thresh=conf_thresh,
  307. cont_thresh=cont_thresh, rp_reject=rp_reject)
  310. def slidingRP_all(spikeTimes, spikeClusters, params=None,
  311. conf_thresh=90, cont_thresh=10, rp_reject=0.0005, n_jobs=1):
  312. """Compute the Sliding RP metric for every cluster in a recording.
  314. :param spikeTimes: array of spike times (s)
  315. :param spikeClusters: array of spike cluster ids that corresponds to spikeTimes
  316. :param params: dict of options passed to slidingRP (e.g. {'recDur': ...})
  317. :param n_jobs: number of parallel processes over clusters. 1 (default) runs
  318. serially; >1 uses that many processes; <0 uses all cores.
  319. (Process-based, like the MATLAB parfor; on Windows the caller
  320. must be under ``if __name__ == '__main__'``.)
  321. :return: dictionary of per-cluster metrics
  322. """
  324. cids = np.unique(spikeClusters)
  325. # Pre-slice spikes per cluster (avoids re-scanning the full arrays per task).
  326. sts = [spikeTimes[spikeClusters == c] for c in cids]
  328. if n_jobs is not None and n_jobs != 1 and len(cids) > 1:
  329. from concurrent.futures import ProcessPoolExecutor
  330. max_workers = None if n_jobs < 0 else n_jobs
  331. args = [(st, params, conf_thresh, cont_thresh, rp_reject) for st in sts]
  332. with ProcessPoolExecutor(max_workers=max_workers) as ex:
  333. results = list(ex.map(_slidingRP_worker, args))
  334. else:
  335. results = [slidingRP(st, params=params, conf_thresh=conf_thresh,
  336. cont_thresh=cont_thresh, rp_reject=rp_reject)
  337. for st in sts]
  339. rpMetrics = {k: [] for k in ('cidx', 'max_confidence', 'min_contamination',
  340. 'rp_min_val', 'n_spikes_below2', 'firing_rate',
  341. 'value', 'value_forced')}
  342. for cid, res in zip(cids, results):
  343. (max_confidence, min_contamination, rp_min_val, n_spikes_below2,
  344. firing_rate, pass_cont_thresh, pass_forced) = res
  345. rpMetrics['cidx'].append(cid)
  346. rpMetrics['max_confidence'].append(max_confidence)
  347. rpMetrics['min_contamination'].append(min_contamination)
  348. rpMetrics['rp_min_val'].append(rp_min_val)
  349. rpMetrics['n_spikes_below2'].append(n_spikes_below2)
  350. rpMetrics['firing_rate'].append(firing_rate)
  351. rpMetrics['value'].append(int(pass_cont_thresh))
  352. rpMetrics['value_forced'].append(int(pass_forced))
  354. return rpMetrics
  358. ## Code from OW
  361. def computeACG(spikeTimes, rpBinSize, nBins):
  362. """Autocorrelogram by histogramming pairwise spike-time differences.
  364. This replicates cortex-lab `histdiff` (used by the MATLAB implementation)
  365. exactly: bin index = floor((t_i - t_j) / rpBinSize), counting each ordered
  366. pair (j < i) whose difference lies in (0, rpBinSize*nBins). Exact
  367. coincidences (diff == 0) are excluded. Implemented with the shifted-train
  368. trick so only nearby pairs are examined.
  370. Parameters
  371. ----------
  372. spikeTimes : numpy.ndarray
  373. Spike times in seconds (need not be sorted).
  374. rpBinSize : float
  375. ACG bin width in seconds (the recording sample period, 1/sampleRate).
  376. nBins : int
  377. Number of ACG bins (the window is rpBinSize * nBins seconds).
  379. Returns
  380. -------
  381. nACG : numpy.ndarray (nBins,)
  382. Integer count of spike pairs with ISI in each bin.
  383. """
  384. st = np.sort(np.asarray(spikeTimes, dtype=np.float64))
  385. max_lag = rpBinSize * nBins
  386. counts = np.zeros(nBins, dtype=np.int64)
  387. n = st.size
  388. shift = 1
  389. # Differences grow with shift, so once none fall within the window we stop.
  390. while shift < n:
  391. dt = st[shift:] - st[:-shift]
  392. m = (dt > 0) & (dt < max_lag)
  393. if not np.any(m):
  394. break
  395. bins = np.floor(dt[m] / rpBinSize).astype(np.int64)
  396. counts += np.bincount(bins, minlength=nBins)[:nBins]
  397. shift += 1
  398. return counts
  401. def computeMatrix(spikeTimes, params):
  402. """Build the [nCont x nRP] confidence matrix for one cluster.
  404. Mirrors matlab/computeMatrix.m. The ACG is computed by computeACG (a
  405. histdiff-equivalent), and expected violations use the Llobet formula with
  406. an explicit recording duration.
  408. Parameters
  409. ----------
  410. spikeTimes : numpy.ndarray
  411. array of spike times (s)
  412. params : dict
  413. - recDur : recording duration (s). Defaults to max(spikeTimes);
  414. recommended to set explicitly.
  415. - binSizeCorr : ACG bin size (s), default 1/sampleRate.
  416. - sampleRate : sample rate (Hz), default 30000.
  417. - cont : vector of contamination levels (%) to test, default 0.5:0.5:35.
  418. - correction : bool, default False. If True, apply the family-wise
  419. multiple-comparisons correction across tau_r (exact Poisson
  420. first-passage; Fig. S3). Slow and over-conservative for short-RP
  421. units; intended for Fig. S3 only.
  422. - rpReject : min tau_r (s) included in the correction, default 0.0005.
  424. Returns
  425. -------
  426. confMatrix : [nCont x nRP] confidence (%) that contamination < cont(i) at rp(j)
  427. cont : tested contamination levels (%)
  428. rp : tested refractory-period durations (s, bin centres)
  429. nACG : ACG counts per bin
  430. firingRate : spike count / recDur (spks/s)
  431. """
  432. sampleRate = params.get('sampleRate', 30000)
  433. rpBinSize = params.get('binSizeCorr', 1 / sampleRate)
  434. recDur = params.get('recDur', None)
  435. if recDur is None:
  436. recDur = np.max(spikeTimes)
  437. # contamination levels (%): 0.5,1,...,35 (70 levels), matching MATLAB
  438. # 0.5:0.5:35 and the manuscript (Python previously stopped at 34.5).
  439. cont = params.get('cont', np.arange(0.5, 35.5, 0.5))
  440. correction = params.get('correction', False)
  441. rp_reject = params.get('rpReject', 0.0005)
  443. rpEdges = np.arange(0, 10 / 1000, rpBinSize) # in s
  444. rp = rpEdges + rpBinSize / 2 # refractory period durations to test (bin centres)
  446. n_spikes = spikeTimes.size
  447. firingRate = n_spikes / recDur
  449. nACG = computeACG(spikeTimes, rpBinSize, rp.size)
  450. obsViol = np.cumsum(nACG[0:rp.size])
  451. refDur = rp + rpBinSize / 2
  453. # Pointwise (nominal) confidence matrix
  454. Nc = n_spikes * (cont[:, np.newaxis] / 100)
  455. Nb = n_spikes * (1 - cont[:, np.newaxis] / 100)
  456. expectedViolMatrix = 2 * refDur[np.newaxis, :] / recDur * Nc * (Nb + (Nc - 1) / 2)
  457. nominalConfMatrix = 100 * (1 - stats.poisson.cdf(obsViol[np.newaxis, :], expectedViolMatrix))
  459. if not correction:
  460. confMatrix = nominalConfMatrix
  461. else:
  462. confMatrix = _fwer_correct(nominalConfMatrix, expectedViolMatrix,
  463. obsViol, rp, rp_reject)
  465. return confMatrix, cont, rp, nACG, firingRate
  468. def _fwer_correct(nominalConfMatrix, expectedViolMatrix, obsViol, rp, rp_reject):
  469. """Family-wise-error-rate correction across tau_r (exact Poisson
  470. first-passage / Markov-chain DP). Port of matlab/computeMatrix.m's
  471. correction branch. Returns the corrected confidence matrix (%)."""
  472. nCont, nRP = nominalConfMatrix.shape
  473. corrected = np.full((nCont, nRP), np.nan)
  474. validIdx = np.where(rp > rp_reject)[0]
  475. if validIdx.size == 0:
  476. corrected[:] = 0
  477. return corrected
  479. for cidx in range(nCont):
  480. expectedViol = expectedViolMatrix[cidx, :]
  481. nomPvals = 1 - nominalConfMatrix[cidx, :] / 100
  482. minPval = np.min(nomPvals[validIdx])
  484. # Shortcut: the corrected confidence cannot exceed (1 - minPval)
  485. if minPval > 0.5:
  486. corrected[cidx, :] = (1 - minPval) * 100
  487. continue
  489. # Per-bin passing boundary: max observed count with pointwise p <= minPval
  490. c = -np.ones(nRP, dtype=np.int64)
  491. expV_valid = expectedViol[validIdx]
  492. obsV_valid = obsViol[validIdx]
  493. max_obs = int(np.max(obsV_valid))
  494. if max_obs == 0:
  495. c[validIdx] = (np.exp(-expV_valid) <= minPval + 1e-10).astype(np.int64) - 1
  496. else:
  497. v_vec = np.arange(max_obs + 1)[:, np.newaxis]
  498. cdf = stats.poisson.cdf(v_vec, expV_valid[np.newaxis, :])
  499. valid_mask = (v_vec <= obsV_valid[np.newaxis, :]) & (cdf <= minPval + 1e-10)
  500. c[validIdx] = valid_mask.sum(axis=0) - 1
  501. max_c = int(np.max(c[validIdx]))
  503. # Markov-chain DP for the first-passage probability
  504. if max_c < 0:
  505. correctedConf = 1.0
  506. else:
  507. P_state = np.zeros(max_c + 1)
  508. P_state[0] = 1.0
  509. P_false = 0.0
  510. lambda_all = np.concatenate(([expectedViol[0]], np.diff(expectedViol)))
  511. for k in range(nRP):
  512. lam = lambda_all[k]
  513. if lam > 0:
  514. pmf = stats.poisson.pmf(np.arange(max_c + 1), lam)
  515. P_state = np.convolve(P_state, pmf)[:max_c + 1]
  516. if c[k] >= 0:
  517. P_false += P_state[:c[k] + 1].sum()
  518. P_state[:c[k] + 1] = 0
  519. correctedConf = 1 - P_false
  521. corrected[cidx, :] = correctedConf * 100
  522. return corrected
  525. def computeViol(obsViol, firingRate, spikeCount, refDur, contaminationProp, recDur):
  526. '''Poisson confidence score for a single (refDur, contamination) hypothesis.
  528. Matches matlab/computeViol.m. Expected violations follow the Llobet et al.
  529. (2022) formulation, in which contaminating spikes produce violations both
  530. with base-neuron spikes and with each other:
  532. Ve = 2 * refDur / recDur * Nc * (Nb + (Nc - 1) / 2)
  534. where Nc = C * N_total and Nb = (1 - C) * N_total.
  536. Parameters
  537. ----------
  538. obsViol : int or array
  539. observed number of violations (cumulative ACG count up to refDur).
  540. firingRate : float
  541. accepted for API parity with the MATLAB signature; not used (recDur
  542. is used directly). Pass None.
  543. spikeCount : int
  544. total spike count of the cluster (N_total).
  545. refDur : float or array
  546. refractory period duration tested, tau_r (seconds).
  547. contaminationProp : float or array
  548. hypothesised contamination as a proportion in [0, 1].
  549. recDur : float
  550. recording duration D (seconds).
  552. Returns
  553. -------
  554. confidenceScore : float or array
  555. probability that true contamination is below contaminationProp given
  556. the observed violations, under a Poisson assumption. 1 - Poisson_CDF.
  557. '''
  558. Nc = spikeCount * contaminationProp
  559. Nb = spikeCount * (1 - contaminationProp)
  560. expectedViol = 2 * refDur / recDur * Nc * (Nb + (Nc - 1) / 2)
  561. confidenceScore = 1 - stats.poisson.cdf(obsViol, expectedViol)
  563. return confidenceScore
  566. def compute_min_contamination(obsViol, spikeCount, refDur, rp, recDur,
  567. conf_thresh=90, rp_reject=0.0005, max_contam=35.0):
  568. """Analytical minimum contamination confirmable at the confidence threshold.
  570. Closed-form equivalent of scanning the contamination grid (cf. Ressmeyer's
  571. analytical method, PR #6), specialised to the exact Llobet Ve. For each
  572. tau_r, the critical Poisson rate at which the confidence equals conf_thresh
  573. is obtained via the Poisson<->chi-squared identity:
  575. lambda_crit = chi2.ppf(conf_thresh/100, 2*(V_o + 1)) / 2
  577. Inverting Ve(C) = lambda_crit (Ve = 2*tau/D * C*N * ((1-C)*N + (C*N-1)/2))
  578. gives the smaller root
  580. C_min(tau) = [ (N - 0.5) - sqrt((N - 0.5)^2 - lambda_crit*D/tau) ] / N.
  582. The unit's minimum confirmable contamination is min_tau C_min(tau) over
  583. tau_r > rp_reject.
  585. Parameters
  586. ----------
  587. obsViol : array cumulative observed violations V_o(tau_r).
  588. spikeCount : int total spikes N.
  589. refDur : array tested tau_r (s) = right bin edge.
  590. rp : array bin-centre tau_r (s), used for the rp_reject mask and
  591. the returned rp_min_val.
  592. recDur : float recording duration D (s).
  593. conf_thresh : float confidence threshold (%).
  594. rp_reject : float exclude tau_r at or below this (s).
  595. max_contam : float cap (%); returns NaN if the minimum exceeds it.
  597. Returns
  598. -------
  599. min_cont : float minimum confirmable contamination (%), or NaN.
  600. rp_min_val : float tau_r (s) at which it is achieved, or NaN.
  601. """
  602. N = spikeCount
  603. if N == 0 or recDur <= 0:
  604. return np.nan, np.nan
  605. gamma = conf_thresh / 100
  606. lam = stats.chi2.ppf(gamma, 2 * (obsViol + 1)) / 2
  607. disc = (N - 0.5) ** 2 - lam * recDur / refDur
  608. Cmin = np.full(refDur.shape, np.nan)
  609. ok = disc >= 0
  610. Cmin[ok] = ((N - 0.5) - np.sqrt(disc[ok])) / N * 100 # as a percentage
  612. testTimes = rp > rp_reject
  613. Cmin_test = Cmin[testTimes]
  614. if not np.any(testTimes) or np.all(np.isnan(Cmin_test)):
  615. return np.nan, np.nan
  616. idx = int(np.nanargmin(Cmin_test))
  617. min_cont = float(Cmin_test[idx])
  618. if min_cont > max_contam:
  619. return np.nan, np.nan
  620. rp_min_val = float(rp[testTimes][idx])
  621. return min_cont, rp_min_val
  625. def plot_acg(ax, acg, timeBins, estimatedIdx=None):
  626. # Convert time in milliseconds for plotting
  627. timeBins = timeBins * 1000
  629. if len(acg) > len(timeBins):
  630. acg = acg[0:len(timeBins)]
  631. # TODO change this so the correct bins are taken for whatever values of timeBins
  632. ax.bar(timeBins, acg, width=np.diff(timeBins)[0], alpha=0.5)
  633. ax.set_xlim(0, timeBins[-1])
  635. ax.spines['right'].set_visible(False)
  636. ax.spines['top'].set_visible(False)
  638. if estimatedIdx is not None:
  639. # Plot RP as black vertical line
  640. ax.plot([timeBins[estimatedIdx], timeBins[estimatedIdx]], [0, max(acg)], 'k-')
  642. ax.set_ylabel('Number of spikes')
  643. ax.set_xlabel('Time (ms)')
  645. def plotSigmoid(ax, acg, timeBins, ySigmoid, estimatedIdx, estimatedRP):
  646. plot_acg(ax, acg, timeBins, estimatedIdx=estimatedIdx)
  647. # Plot on top of ACG the sigmoid
  648. ax.plot(timeBins[0:len(ySigmoid)] * 1000, ySigmoid, 'b')
  649. ax.plot(timeBins[estimatedIdx] * 1000, ySigmoid[estimatedIdx], 'rx')
  650. ax.set_title('Estimated RP:%.2f ms' % estimatedRP)
  653. def plotSlidingRP(spikeTimes, params=None, plotXs=None, inputAxes=None,
  654. plotExtraContours=False):
  655. """Visualise the Sliding RP result for one cluster (Python twin of
  656. matlab/plotSlidingRP.m). Produces three panels:
  657. 1. ACG (0-5 ms),
  658. 2. confidence matrix (contamination x tau_r) with the 90% iso-contour and
  659. the minimum-contamination time,
  660. 3. confidence trace at the contamination threshold (10%), optionally also
  661. at 7.5% and 15% (plotExtraContours).
  663. Parameters
  664. ----------
  665. spikeTimes : array of spike times (s).
  666. params : dict, optional. Recognised keys: sampleRate, binSizeCorr, recDur,
  667. contaminationThresh (default 10), confidenceThresh (default 90),
  668. savefig (bool), figpath (str, no extension).
  669. plotXs : [taus_ms, colors] optional. Vertical markers at the given tau_r
  670. values (ms) in the matching colors, drawn on panels 1 and 3.
  671. inputAxes : (fig, axs) optional, where axs has at least 3 axes; otherwise a
  672. new 1x3 figure is created.
  673. plotExtraContours : bool. Also draw the 7.5% and 15% confidence traces.
  675. Returns the (fig, axs) used.
  676. """
  677. import matplotlib.pyplot as plt
  679. params = dict(params) if params else {}
  680. contThresh = params.get('contaminationThresh', 10)
  681. confThresh = params.get('confidenceThresh', 90)
  683. confMatrix, cont, rp, nACG, firingRate = computeMatrix(np.asarray(spikeTimes), params)
  684. _, _, rp_min_val, _, _, _, _ = slidingRP(
  685. spikeTimes, params=params, conf_thresh=confThresh, cont_thresh=contThresh)
  687. if inputAxes is not None:
  688. fig, axs = inputAxes
  689. else:
  690. fig, axs = plt.subplots(1, 3, figsize=(12, 4))
  692. rp_ms = rp * 1000
  694. # --- Panel 1: ACG ---
  695. ax = axs[0]
  696. ax.bar(rp_ms, nACG[0:len(rp)], width=np.diff(rp_ms)[0], color='k')
  697. ax.set_xlim([0, 5])
  698. ax.set_xlabel('Time from spike (ms)')
  699. ax.set_ylabel('ACG count (spks)')
  700. ax.fill(np.array([0, 1, 1, 0]) * 0.5, np.array([0, 0, 1, 1]) * ax.get_ylim()[1],
  701. 'k', alpha=0.2)
  702. if plotXs is not None:
  703. for tau_ms, c in zip(plotXs[0], plotXs[1]):
  704. ax.axvline(tau_ms, color=c, linewidth=1)
  705. ax.spines['right'].set_visible(False)
  706. ax.spines['top'].set_visible(False)
  708. # --- Panel 2: confidence matrix ---
  709. ax = axs[1]
  710. im = ax.imshow(confMatrix, extent=[rp_ms[0], rp_ms[-1], cont[0], cont[-1]],
  711. aspect='auto', vmin=0, vmax=100, origin='lower')
  712. cbar = fig.colorbar(im, ax=ax)
  713. cbar.set_label('Confidence (%)')
  714. ax.plot([rp_ms[0], rp_ms[-1]], [contThresh, contThresh], 'r', linewidth=1)
  715. if not np.isnan(rp_min_val):
  716. ax.plot([rp_min_val * 1000] * 2, [cont[0], cont[-1]], 'r', linewidth=1)
  717. # 90% iso-contour
  718. z = np.vstack([np.zeros((1, confMatrix.shape[1])), confMatrix])
  719. ii = np.argmax(z > confThresh, axis=0).astype(float)
  720. ii[ii == 0] = np.nan
  721. contContour = np.full(ii.shape, np.nan)
  722. good = ~np.isnan(ii)
  723. contContour[good] = cont[(ii[good] - 1).astype(int)]
  724. ax.plot(rp_ms, contContour, 'r', linewidth=2)
  725. ax.set_xlim([0, 5])
  726. ax.set_xlabel('Time from spike (ms)')
  727. ax.set_ylabel('Contamination (%)')
  728. ax.invert_yaxis()
  730. # --- Panel 3: confidence trace(s) ---
  731. ax = axs[2]
  732. def _trace(level, **kw):
  733. idx = np.argmin(np.abs(cont - level))
  734. ax.plot(rp_ms, confMatrix[idx, :], **kw)
  735. if plotExtraContours:
  736. _trace(7.5, color='lightcoral', linewidth=1, label='7.5%')
  737. _trace(15, color='darkred', linewidth=1, label='15%')
  738. _trace(contThresh, color='r', linewidth=2, label='%g%%' % contThresh)
  739. ax.plot([0, 5], [confThresh, confThresh], 'k', linewidth=1)
  740. if plotXs is not None:
  741. for tau_ms, c in zip(plotXs[0], plotXs[1]):
  742. ax.axvline(tau_ms, color=c, linewidth=1)
  743. ax.set_xlim([0, 5])
  744. ax.set_ylim([0, 100])
  745. ax.set_xlabel('Time from spike (ms)')
  746. ax.set_ylabel('Confidence of <=%g%% contamination (%%)' % contThresh)
  747. ax.legend(frameon=False, fontsize=8)
  748. ax.spines['right'].set_visible(False)
  749. ax.spines['top'].set_visible(False)
  751. fig.tight_layout()
  752. if params.get('savefig', False) and 'figpath' in params:
  753. fig.savefig(params['figpath'] + '.svg', dpi=300)
  754. fig.savefig(params['figpath'] + '.png', dpi=300)
  755. return fig, axs
  757. # return acg
  758. #
  759. # # helper functions
  760. #
  761. #
  762. # def find_nearest(array, value):
  763. # array = np.asarray(array)
  764. # subtracted = (array - value)
  765. # valid_idx = np.where(subtracted >= 0)[0]
  766. # if len(valid_idx) > 0:
  767. #
  768. # out = valid_idx[subtracted[valid_idx].argmin()]
  769. #
  770. # else:
  771. # out = np.nan
  772. #
  773. # return out
  774. #
  775. #

metrics.py at commit bb841e5, no license · at the source

Overview

Authors: Noam Roth1, Gaelle Chapuis2, Olivier Winter2, International Brain Laboratory2, Ryan A. Ressmeyer1,3, Luke M. Bun1,3, Ryan A. Canfield1,3, Gregory D. Horwitz1, Nicholas A. Steinmetz1
ORCID iDs: Noam Roth, Gaelle Chapuis, Olivier Winter, Ryan A. Ressmeyer, Luke M. Bun, Ryan A. Canfield, Gregory D. Horwitz, Nicholas A. Steinmetz
  1. Department of Neurobiology & Biophysics, University of Washington, Seattle, WA, USA
  2. International Brain Laboratory
  3. Department of Bioengineering, University of Washington, Seattle, WA, USA
Institutions: University of Washington (United States); International Brain Laboratory (Portugal)
Journal: bioRxiv (preprint)
Dates: published online 9 March 2026
Type: Preprint · Language: English
License: CC BY-NC
Identifiers: DOI 10.64898/2026.03.06.710130 · OpenAlex W7135018293
Open access: green, a free copy (OpenAlex)
Status: code verified
Categories: methods / tools (subfield)
Methods: Statistics, Preprocessing, Connectivity, fMRI & imaging, Single-unit activity, calcium imaging, Physiology & signal measures
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Wellcome Trust (216324/Z/19/Z)
Citations: not cited yet (Europe PMC); 27 references in the paper

Abstract

The increasing size of electrophysiological datasets has heightened the need for quality metrics that automatically reject neurons whose activity was recorded with low sensitivity or specificity. One key approach estimates artifactual contamination by assuming that each neuron has a refractory period (RP), a brief time interval following each action potential when further activity cannot occur. However, existing methods cannot be applied without prior knowledge of the neurons’ RP durations, limiting their usefulness in datasets that include neurons from brain regions or species in which RP durations have not been systematically characterized. Here, we find that neurons in some brain regions (thalamus) and species (macaque) have shorter RP durations than commonly assumed, and we introduce a new metric, the Sliding Refractory Period metric, which is robust to variation in a neuron’s RP duration without tuning. We validate the method using simulations, demonstrating that it improves acceptance of uncontaminated spike trains with short or long RP durations while still rejecting contaminated ones. Moreover, by incorporating Poisson statistics into the calculation, the method also improves on prior work by allowing the user to approximately control the false acceptance rate. Our new metric improves quantification of contamination in electrophysiological recordings and enables application of a single tuning-free quality metric to data recorded from diverse brain regions and species.

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

Repository

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

SteinmetzLab/slidingRefractory

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: bb841e55defb9cb95242a8177d50de1916c2a2e2, 16 July 2026
Languages: MATLAB (24), Python (15)
Size: 62 files, 39 scripts
Software Heritage: not archived
Found in: “Code and Data Availability”
Holds: README, CITATION.cff, environment (pyproject.toml), tests, continuous integration
Not found: license file, documentation
Tools: NumPy (12 files), Matplotlib (6 files), pandas (4 files), SciPy (4 files), Signal Processing Toolbox (2 files), Optimization Toolbox (1 file), Parallel Computing Toolbox (1 file), Statistics and Machine Learning Toolbox (1 file), Phy (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
40 files
At the source: github.com/SteinmetzLab/slidingRefractory/

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;
  • 39 scripts, each with its path and the digest of its content;
  • 16 matches 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.

Code and Data Availability

Code for the metric and for the simulations presented is available at https://github.com/SteinmetzLab/slidingRefractory/.

Reproduced under the paper's license (CC BY-NC), 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, 30 September 2026: the first record

Recorded: type, language, journal, dates, 9 authors, 1 funder, 24 references.

Cite

This paper

Roth, N., Chapuis, G., Winter, O., International Brain Laboratory, Ressmeyer, R. A., Bun, L. M., Canfield, R. A., Horwitz, G. D., & Steinmetz, N. A. (2026). A flexible quality metric for electrophysiological recordings across brain regions and species. bioRxiv (preprint). https://doi.org/10.64898/2026.03.06.710130

BibTeX

@article{roth2026flexible,
author = {Roth, Noam and Chapuis, Gaelle and Winter, Olivier and {International Brain Laboratory} and Ressmeyer, Ryan A. and Bun, Luke M. and Canfield, Ryan A. and Horwitz, Gregory D. and Steinmetz, Nicholas A.},
title = {{A flexible quality metric for electrophysiological recordings across brain regions and species}},
journal = {bioRxiv (preprint)},
year = {2026},
month = mar,
publisher = {bioRxiv},
issn = {2692-8205},
doi = {10.64898/2026.03.06.710130},
url = {https://doi.org/10.64898/2026.03.06.710130}
}

RIS

TY - JOUR
AU - Roth, Noam
AU - Chapuis, Gaelle
AU - Winter, Olivier
AU - International Brain Laboratory
AU - Ressmeyer, Ryan A.
AU - Bun, Luke M.
AU - Canfield, Ryan A.
AU - Horwitz, Gregory D.
AU - Steinmetz, Nicholas A.
TI - A flexible quality metric for electrophysiological recordings across brain regions and species
T2 - bioRxiv (preprint)
J2 - bioRxiv
PY - 2026
DA - 2026/03/09
SN - 2692-8205
PB - bioRxiv
DO - 10.64898/2026.03.06.710130
UR - https://doi.org/10.64898/2026.03.06.710130
LA - en
ER -

CSL-JSON

{
"id": "10.64898/2026.03.06.710130",
"type": "article",
"title": "A flexible quality metric for electrophysiological recordings across brain regions and species",
"container-title": "bioRxiv (preprint)",
"author": [
{
"family": "Roth",
"given": "Noam"
},
{
"family": "Chapuis",
"given": "Gaelle"
},
{
"family": "Winter",
"given": "Olivier"
},
{
"literal": "International Brain Laboratory"
},
{
"family": "Ressmeyer",
"given": "Ryan A."
},
{
"family": "Bun",
"given": "Luke M."
},
{
"family": "Canfield",
"given": "Ryan A."
},
{
"family": "Horwitz",
"given": "Gregory D."
},
{
"family": "Steinmetz",
"given": "Nicholas A."
}
],
"container-title-short": "bioRxiv",
"DOI": "10.64898/2026.03.06.710130",
"ISSN": "2692-8205",
"publisher": "bioRxiv",
"URL": "https://doi.org/10.64898/2026.03.06.710130",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
9
]
]
}
}

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.1016/j.patter.2026.101590 [code]
Density-based longitudinal neuron tracking in high-density electrophysiological recordings.
Journal: Patterns (New York, N.Y.)
In common: Phy, Optimization Toolbox, Parallel Computing Toolbox, 7 other tools, 3 references
[2] doi:10.1038/s41593-026-02357-2 [code]
Experience reorganizes content-specific memory traces in macaques.
Journal: Nature neuroscience
In common: Phy, Optimization Toolbox, Parallel Computing Toolbox, 6 other tools, 1 reference
[3] doi:10.1038/s41467-026-75347-4 [code]
Sleep reveals dynamics integrating and segregating movement and stimulus representations in V1.
Journal: Nature communications
In common: Phy, Optimization Toolbox, Parallel Computing Toolbox, 6 other tools, 1 reference
[4] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: Phy, Optimization Toolbox, Signal Processing Toolbox, 6 other tools, 1 reference
[5] doi:10.1038/s41467-026-71331-0 [code]
A multimodal approach for visualizing and identifying electrophysiological cell types in vivo.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, 3 references, author Nicholas A. Steinmetz
[6] doi:10.1038/s41592-026-03076-z [code]
Neuropixels Opto: combining high-resolution electrophysiology and optogenetics.
Journal: Nature methods
In common: pandas, SciPy, Matplotlib, 1 other tool, 3 references, author Nicholas A. Steinmetz
[7] doi:10.1038/s41467-026-71725-0 [code]
Interactions across hemispheres in prefrontal cortex reflect global cognitive processing.
Journal: Nature communications
In common: Optimization Toolbox, Parallel Computing Toolbox, Signal Processing Toolbox, 5 other tools, 2 references
[8] doi:10.1523/jneurosci.2001-25.2026 [code]
Dynamics of Dentate Gyrus Place Cells and Dentate Spikes during Spatial and Nonspatial Changes in Environments.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: Phy, Optimization Toolbox, Signal Processing Toolbox, 5 other tools, 1 reference
[9] doi:10.1038/s41467-026-71664-w [code]
Dorsal prefrontal cortex drives perseverative behavior in mice.
Journal: Nature communications
In common: Optimization Toolbox, Parallel Computing Toolbox, Signal Processing Toolbox, 2 other tools, 4 references
[10] doi:10.1126/sciadv.aeh7220 [code]
Central complex representations of self-movement are sufficient to compute wind direction in flight.
Journal: Science advances
In common: Optimization Toolbox, Parallel Computing Toolbox, Signal Processing Toolbox, 6 other tools

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.