OSCR

Hippocampal skill memory expansion drives online performance dynamics during skill learning.

Code ↔ Paper

2 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 2 matches
  1. [1] § STAR★METHODS › METHOD DETAILS › Oscillatory bursting analysis ↔ spectralevents.py, lines 293–399 · score 0.79 · spectralEvents, full width, median power, Event duration, FWHM, detection
  2. [2] § STAR★METHODS › METHOD DETAILS › Oscillatory bursting analysis ↔ spectralevents_find.m, lines 282–404 · score 0.79 · spectralEvents, full width, median power, Event duration, FWHM, overlap

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 · 601 lines · 23 KB · BSD-3-Clause · 1 match

  1. '''Spectral Event Analysis functions'''
  3. # Authors: Tim Bardouille <[email hidden]>
  4. # Ryan Thorpe <[email hidden]>
  6. import numpy as np
  7. import scipy.signal as signal
  8. import scipy.ndimage as ndimage
  9. import scipy.ndimage.filters as filters
  10. import matplotlib.pyplot as plt
  13. def tfr(timeseries, freqs, samp_freq, width=7.):
  14. '''Calculate the time-frequency response by convolving with Morlet wavelets
  16. Parameters
  17. ----------
  18. timeseries : array-like, shape ([n_epochs,] n_times)
  19. The timeseries signal, one for each epoch of data.
  20. freqs : array-like, shape (n_freqs,)
  21. Frequency domain of the TFR (Hz).
  22. samp_freq : float
  23. Sampling frequency (Hz).
  24. width : int
  25. Number of cycles in each Morlet wavelet (>5 advisable).
  27. Returns
  28. -------
  29. tfr : array, shape (n_epochs, n_freqs, n_times)
  30. A collection of time-frequency responses (TFRs), one for each epoch.
  32. Notes
  33. -----
  34. Adapted from Ole Jensen's traces2tfr in the 4Dtools toolbox.
  35. '''
  37. ts = np.atleast_2d(timeseries)
  38. n_epochs = ts.shape[0]
  39. n_samps = ts.shape[1]
  40. n_freqs = len(freqs)
  42. # Validate freqs input
  43. f_nyquist = samp_freq / 2 # Nyquist frequency
  44. dt = 1 / samp_freq # Sampling time interval
  45. min_freq = 1 / (n_samps * dt) # Minimum resolvable frequency
  47. if freqs[0] < min_freq:
  48. raise ValueError('Frequency vector includes values outside the '
  49. 'resolvable/alias-free range.')
  50. elif freqs[-1] > f_nyquist:
  51. raise ValueError('Frequency vector includes values outside the '
  52. 'resolvable/alias-free range.')
  53. elif np.abs(freqs[1] - freqs[0]) < min_freq:
  54. raise ValueError('Frequency vector includes values outside the '
  55. 'resolvable/alias-free range.')
  57. tfr = np.zeros((n_epochs, n_freqs, n_samps))
  59. # Trial Loop
  60. for trial_idx in np.arange(n_epochs):
  61. ts_detrended = signal.detrend(ts[trial_idx, :])
  62. # Frequency loop
  63. for freq_idx in np.arange(n_freqs):
  64. tfr[trial_idx, freq_idx, :] = _energyvec(freqs[freq_idx],
  65. ts_detrended, samp_freq,
  66. width)
  68. return tfr
  71. def _get_power_thresholds(tfr, FOM_threshold=6.):
  72. '''Get the power threshold for each frequency band of a TFR'''
  74. med_powers = np.median(tfr, axis=1)
  75. return med_powers * FOM_threshold, med_powers
  78. def tfr_normalize(tfr):
  79. '''Normalize the power in each frequency band of a TFR
  81. Parameters
  82. ----------
  83. tfr : array, shape ([n_epochs,] n_freqs, n_times)
  84. The time-frequency response (TFR) to be normalized.
  86. Returns
  87. -------
  88. tfr_norm : array
  89. The normalized TFR calculated by dividing the power values in each
  90. frequency bin by the median power across all trials and time samples.
  91. '''
  93. if len(tfr.shape) == 3:
  94. n_epochs, _, n_times = tfr.shape
  95. med_powers = np.median(tfr, axis=(0, 2))
  96. med_powers_tiled = np.tile(med_powers, (n_epochs, n_times, 1))
  97. med_powers = np.transpose(med_powers_tiled, axes=(0, 2, 1))
  99. elif len(tfr.shape) == 2:
  100. _, n_times = tfr.shape
  101. med_powers = np.median(tfr, axis=1)
  102. med_powers_tiled = np.tile(med_powers, (n_times, 1))
  103. med_powers = np.transpose(med_powers_tiled, axes=(1, 0))
  105. else:
  106. raise ValueError(f'TFR must be an array of at least 2 dimensions. Got '
  107. f'{tfr.shape}.')
  109. return tfr / med_powers
  112. def find_events(tfr, times, freqs, event_band, thresholds=None,
  113. threshold_FOM=6.):
  114. '''Locate spectral events in a time-frequency response.
  116. Parameters
  117. ----------
  118. tfr : array, shape ([n_epochs,] n_freqs, n_times)
  119. The time-frequency response (TFR) in which to search for high power
  120. spectral events.
  121. times : array-like, shape (n_times,)
  122. Time domain of the TFR in seconds.
  123. freqs : array-like, shape (n_freqs,)
  124. Frequency domain of the TFR in Hertz.
  125. event_band : list
  126. Lower and upper bounds (inclusive, Hz) of the frequency band-of-
  127. interest in which to search for spectral events.
  128. thresholds : None | array-like, shape (n_freqs,)
  129. If not None, these frequency-specific threshold values (in units of
  130. spectral power) will be used to identify suprathreshold spectral
  131. events.
  132. threshold_FOM : float | int
  133. Factor-of-the-median threshold (a.u.) with which to identify
  134. suprathreshold spectral events (default: 6). This threshold value is
  135. applied across frequencies of the TFR when thresholds=None. See Shin et
  136. al. eLife 2017 for more details concerning this value.
  138. Returns
  139. -------
  140. events : list of list of dict
  141. A nested list with n_epochs elements in the outer list and n_events in
  142. the inner list. Each inner element comprises an event that is
  143. characterized by dictionary items including it's location in time,
  144. location in frequency, duration, and frequency span, and more.
  146. Notes
  147. -----
  148. This version only supports find-method #1 at the moment. As outlined in
  149. Shin et al. eLife (2017), it isolates spectral events by first retrieving
  150. all local maxima in un-normalized TFR using imregionalmax, then selecting
  151. suprathreshold peaks within the frequency band of interest. This method
  152. allows for multiple, overlapping events to occur in a given suprathreshold
  153. region and does not guarantee that the presence of within-band,
  154. suprathreshold activity in any given trial will render an event.
  155. '''
  157. # ensure tfr has 3 dimensions (epochs x freqs x time)
  158. if len(tfr.shape) < 3:
  159. tfr = tfr[np.newaxis, ...]
  160. n_epochs = tfr.shape[0]
  161. n_freqs = tfr.shape[1]
  162. n_times = tfr.shape[2]
  164. # some time steps might be slightly different due to rounding error, etc.
  165. samp_freq = 1 / np.unique(np.diff(times).round(10))
  166. if len(samp_freq) > 1:
  167. raise ValueError('Sampling rate is not consistent across time '
  168. 'samples.')
  169. samp_freq = samp_freq[0]
  171. # concatenate trials together to make one big spectrogram, then find thresh
  172. tfr_permute = np.transpose(tfr, [1, 2, 0]) # freq x time x trial
  173. tfr_concat_trials = np.reshape(tfr_permute, (n_freqs, n_times * n_epochs))
  174. thresholds, med_powers = _get_power_thresholds(tfr_concat_trials,
  175. FOM_threshold=threshold_FOM)
  177. # Validate consistency of parameter dimensions
  178. if n_freqs != len(freqs):
  179. raise ValueError('Mismatch in frequency dimensions!')
  180. if n_times != len(times):
  181. raise ValueError('Mismatch in time dimensions!')
  183. # Find spectral events using appropriate method
  184. # Implementing find_method=1 for now
  185. events = _find_localmax_method_1(tfr, freqs, times, event_band,
  186. thresholds, med_powers, samp_freq)
  188. return events
  191. def _energyvec(f, s, Fs, width=7.):
  192. '''Spectral energy as a function of frequency using Morlet wavelets.'''
  194. dt = 1 / Fs
  195. sf = f / width
  196. st = 1 / (2 * np.pi * sf)
  198. t_single_side = np.arange(0., 3.5 * st, dt)
  199. # mirror about 0; always an odd # of elements
  200. times = np.r_[-t_single_side[-1:0:-1], t_single_side]
  201. wavelet = _morlet(f, times, width)
  203. fourier_coeffs = np.convolve(s, wavelet)
  204. spec_energy = 2 * (dt * np.abs(fourier_coeffs)) ** 2
  205. lower_idx = int(len(wavelet) // 2)
  206. upper_idx = int(len(spec_energy) - (len(wavelet) // 2))
  208. return spec_energy[lower_idx:upper_idx]
  211. def _morlet(f, t, width):
  212. '''
  213. Morlet's wavelet for frequency f and time t. The wavelet will be normalized
  214. so the total energy is 1. width defines the ``width'' of the wavelet. A
  215. value >= 5 is suggested.
  217. Ref: Tallon-Baudry et al., J. Neurosci. 15, 722-734 (1997)
  218. '''
  220. sf = f / width
  221. st = 1 / (2 * np.pi * sf)
  222. A = 1 / (st * np.sqrt(2 * np.pi))
  223. y = A * np.exp(-t ** 2 / (2 * st ** 2)) * np.exp(1j * 2 * np.pi * f * t)
  225. return y
  228. def _fwhm_lower_upper_bound1(vec, peakInd, peakValue):
  229. '''
  230. Function to find the lower and upper indices within which the vector is
  231. less than the FWHM with some rather complicated boundary rules
  232. (Shin, eLife, 2017).
  233. '''
  235. halfMax = peakValue/2
  237. # Extract data before the peak only (data should be rising at the end of
  238. # the new array)
  239. vec1 = vec[0:peakInd]
  240. # Find indices less than half the max
  241. vec1_underThreshold = np.where(vec1 < halfMax)[0]
  242. if len(vec1_underThreshold) == 0:
  243. # There are no indices less than half the max, so we have to estimate
  244. # the lower edge
  245. estimateLowerEdge = True
  246. else:
  247. # There are indices less than half the max, take the last one under
  248. # halfMax as the lower edge
  249. estimateLowerEdge = False
  250. lowerEdgeIndex = vec1_underThreshold[-1]
  252. # Extract data following the peak only (data should be falling at the start
  253. # of the new array)
  254. vec2 = vec[peakInd:]
  255. # Find indices less than half the max
  256. vec2_underThreshold = np.where(vec2 < halfMax)[0]
  257. if len(vec2_underThreshold) == 0:
  258. # There are no indices less than half the max, so we have to estimate
  259. # the upper edge
  260. estimateUpperEdge = True
  261. else:
  262. # There are indices less than half the max, take the first one under
  263. # halfMax as the upper edge
  264. estimateUpperEdge = False
  265. upperEdgeIndex = vec2_underThreshold[0] + len(vec1)
  267. if not estimateLowerEdge:
  268. if not estimateUpperEdge:
  269. # FWHM fits in the range, so pick off the edges of the FWHM
  270. lowerInd = lowerEdgeIndex
  271. upperInd = upperEdgeIndex
  272. FWHM = upperInd - lowerInd
  273. if estimateUpperEdge:
  274. # FWHM fits in on the low end, but hits the edge on the high end
  275. lowerInd = lowerEdgeIndex
  276. upperInd = len(vec)-1
  277. FWHM = 2 * (peakInd - lowerInd + 1)
  278. else:
  279. if not estimateUpperEdge:
  280. # FWHM hits the edge on the low end, but fits on the high end
  281. lowerInd = 0
  282. upperInd = upperEdgeIndex
  283. FWHM = 2 * (upperInd - peakInd + 1)
  284. if estimateUpperEdge:
  285. # FWHM hits the edge on the low end and the high end
  286. lowerInd = 0
  287. upperInd = len(vec)-1
  288. FWHM = 2*len(vec)
  290. return lowerInd, upperInd, FWHM
  293. def _find_localmax_method_1(tfr, freqs, times, event_band,
  294. eventThresholdByFrequency,
  295. medianPower, Fs):
  296. '''
  297. 1st event-finding method (primary event detection method in Shin et
  298. al. eLife 2017): Find spectral events by first retrieving all local
  299. maxima in un-normalized TFR using imregionalmax, then selecting
  300. suprathreshold peaks within the frequency band of interest. This
  301. method allows for multiple, overlapping events to occur in a given
  302. suprathreshold region and does not guarantee the presence of
  303. within-band, suprathreshold activity in any given trial will render
  304. an event.
  306. spectralEvents: 12 column matrix for storing local max event metrics:
  307. hit/miss, maxima frequency,
  308. lowerbound frequency, upperbound frequency,
  309. frequency span, maxima timing, event onset timing,
  310. event offset timing, event duration, maxima power,
  311. maxima/median power
  312. '''
  314. n_epochs = tfr.shape[0]
  316. all_epochs_events = list()
  318. # Retrieve all local maxima in tfr using python equivalent of imregionalmax
  319. for trial_idx in range(n_epochs):
  321. epoch_events = list()
  323. # Get tfr data for this trial [frequency x time]
  324. thistfr = tfr[trial_idx, :, :]
  326. # Find local maxima in the tfr data
  327. data = thistfr
  328. # Find maximum amoung adjacent pixels (3x3 footprint) for each pixel
  329. data_max = filters.maximum_filter(data, size=(3, 3))
  330. maxima = (data == data_max)
  331. data_min = filters.minimum_filter(data, size=(3, 3))
  332. # Rule out pixels with footprints that have flatlined
  333. maxima[data_max == data_min] = False
  334. labeled, num_objects = ndimage.label(maxima)
  335. yx = ndimage.center_of_mass(data, labels=labeled,
  336. index=range(1, num_objects + 1))
  338. peakF = list()
  339. peakT = list()
  340. peakPower = list()
  341. for f_idx, t_idx in yx:
  342. f_idx = int(round(f_idx))
  343. t_idx = int(round(t_idx))
  344. event_freq = freqs[f_idx]
  345. if (event_freq >= event_band[0] and event_freq <= event_band[1] and
  346. thistfr[f_idx, t_idx] > eventThresholdByFrequency[f_idx]):
  347. peakF.append(f_idx)
  348. peakT.append(t_idx)
  349. peakPower.append(thistfr[f_idx, t_idx])
  351. numPeaks = len(peakF)
  353. # Find local maxima lowerbound, upperbound, and full width at half max
  354. # for both frequency and time
  355. for lmi in range(numPeaks):
  356. thisPeakF = peakF[lmi]
  357. thisPeakT = peakT[lmi]
  358. thisPeakPower = peakPower[lmi]
  360. # Indices of tfr frequencies < half max power at the time of a
  361. # given local peak
  362. tfrFrequencies = thistfr[:, thisPeakT]
  363. lowerInd, upperInd, FWHM = _fwhm_lower_upper_bound1(tfrFrequencies,
  364. thisPeakF,
  365. thisPeakPower)
  366. lowerEdgeFreq = freqs[lowerInd]
  367. upperEdgeFreq = freqs[upperInd]
  368. FWHMFreq = FWHM * (freqs[1] - freqs[0])
  370. # Indices of tfr times < half max power at the frequency of a given
  371. # local peak
  372. tfrTimes = thistfr[thisPeakF, :]
  373. lowerInd, upperInd, FWHM = _fwhm_lower_upper_bound1(tfrTimes,
  374. thisPeakT,
  375. thisPeakPower)
  376. lowerEdgeTime = times[lowerInd]
  377. upperEdgeTime = times[upperInd]
  378. FWHMTime = FWHM / Fs
  380. # Put peak characteristics to a dictionary
  381. peakParameters = {
  382. 'Peak Frequency': freqs[thisPeakF],
  383. 'Lower Frequency Bound': lowerEdgeFreq,
  384. 'Upper Frequency Bound': upperEdgeFreq,
  385. 'Frequency Span': FWHMFreq,
  386. 'Peak Time': times[thisPeakT],
  387. 'Event Onset Time': lowerEdgeTime,
  388. 'Event Offset Time': upperEdgeTime,
  389. 'Event Duration': FWHMTime,
  390. 'Peak Power': thisPeakPower,
  391. 'Normalized Peak Power': thisPeakPower / medianPower[thisPeakF]
  392. }
  394. # Build a list of dictionaries
  395. epoch_events.append(peakParameters)
  397. all_epochs_events.append(epoch_events)
  399. return all_epochs_events
  402. def plot_events(tfr, times, freqs, event_band, spec_events=None,
  403. timeseries=None, ax=None, vlim=None, ylim_ts=None, label=None):
  404. '''Plot a single TFR spectrogram overlayed with spectral events.
  406. Parameters
  407. ----------
  408. tfr : array, shape (n_freqs, n_times)
  409. A time-frequency response (TFR).
  410. times : array-like, shape (n_times,)
  411. Time domain of the TFR in seconds.
  412. freqs : array-like, shape (n_freqs,)
  413. Frequency domain of the TFR in Hertz.
  414. event_band : list
  415. Lower and upper bounds (inclusive, Hz) of the frequency band-of-
  416. interest in which spectral events were detected.
  417. spec_events : list of dict
  418. A single-level list where each element comprises an event that is
  419. characterized by dictionary items including it's location in time,
  420. location in frequency, duration, and frequency span, and more. If None,
  421. no events are plotted.
  422. timeseries : array, shape (n_times,) | None
  423. The timeseries associated with the TFR.
  424. ax : None | Matplotlib.axes.Axes instance
  425. The Axes instance with wich to plot the spectrogram on. If None, a new
  426. Axes object is created.
  427. vlim : list, shape (2,) | None
  428. If not None, sets the colorbar lower and upper bounds for the
  429. spectrogram across all axes of the returned figure. If None, each
  430. spectrogram (including example epochs) sets its vlim independently.
  431. ylim_ts : list | None
  432. The lower and upper limits of the timeseries (if applicable) overlayed
  433. on the spectrogram. If None, these values are automatically assigned.
  434. label : str | None
  435. A label to overlay in the upper left corner of this plot.
  437. Return
  438. ------
  439. fig : Matplotlib.figure.Figure instance
  440. The Figure instance containing the spectrogram.
  441. '''
  443. if tfr.shape != (len(freqs), len(times)):
  444. raise ValueError(f'tfr must be an array of shape (n_freqs, n_times), '
  445. f'got tfr: {tfr.shape}, freqs: ({len(freqs)},), '
  446. f'times: ({len(times)},)')
  448. if vlim is None:
  449. vlim = [None, None]
  451. # convert to numpy array if not already
  452. freqs = np.array(freqs)
  454. # frequencies within the band of interest
  455. # band_mask = np.logical_and(freqs >= event_band[0],
  456. # freqs <= event_band[1])
  458. if ax is None:
  459. fig, ax = plt.subplots(1, 1)
  460. else:
  461. fig = ax.get_figure()
  463. # plot tfr
  464. im = ax.pcolormesh(times, freqs, tfr, cmap='jet', vmin=vlim[0],
  465. vmax=vlim[1], shading='nearest')
  466. fig.colorbar(im, ax=ax)
  467. ax.axhline(y=event_band[0], c='w', linewidth=2., linestyle=':', alpha=.7)
  468. ax.axhline(y=event_band[1], c='w', linewidth=2., linestyle=':', alpha=.7)
  469. ax.set_yticks([freqs[0], event_band[0], event_band[1], freqs[-1]])
  470. ax.set_xlim(times[0], times[-1])
  472. # overlay with timecourse
  473. # twin axis needs to be created and yticks set regardless of if a
  474. # timeseries is provided to ensure consistency with seaborn tick formatting
  475. ax_twin = ax.twinx()
  476. ax_twin.set_yticks([])
  477. if timeseries is not None:
  478. ax_twin.plot(times, timeseries, 'w', linewidth=1., alpha=0.8)
  480. if ylim_ts is not None:
  481. ax_twin.set_ylim(ylim_ts[0], ylim_ts[1])
  483. # plot event locations
  484. if spec_events is not None:
  485. event_times = [event['Peak Time'] for event in spec_events]
  486. event_freqs = [event['Peak Frequency'] for event in spec_events]
  487. # reverse sigmoid: make scatter markers more transparent when there are
  488. # more of them
  489. alpha = lambda x : (0.6 + 0.4 * np.exp(-0.2 * (x - 30)) # noqa
  490. / (1 + np.exp(-0.2 * (x - 30)))) # noqa
  491. ax.scatter(event_times, event_freqs, s=20, c='w', marker='x',
  492. alpha=alpha(len(spec_events)))
  494. # add label on top of spectrogram plot
  495. if label is not None:
  496. ax.annotate(label, xy=(0.01, 0.98), xycoords='axes fraction',
  497. va='top', ha='left', color='w', size=10, fontweight='bold')
  499. return fig
  502. def plot_avg_spectrogram(tfr, times, freqs, event_band, spec_events=None,
  503. timeseries=None, example_epochs=None, vlim=None,
  504. show_events=False):
  505. '''Plot the average spectrogram over all TFR epochs/trials.
  507. Parameters
  508. ----------
  509. tfr : array, shape (n_epochs, n_freqs, n_times)
  510. A stack of time-frequency response (TFR) epochs/trials.
  511. times : array-like, shape (n_times,)
  512. Time domain of each TFR in seconds.
  513. freqs : array-like, shape (n_freqs,)
  514. Frequency domain of each TFR in Hertz.
  515. event_band : list
  516. Lower and upper bounds (inclusive, Hz) of the frequency band-of-
  517. interest in which spectral events were detected.
  518. spec_events : list of list of dict
  519. A nested list with n_epochs elements in the outer list and n_events in
  520. the inner list. Each inner element comprises an event that is
  521. characterized by dictionary items including it's location in time,
  522. location in frequency, duration, and frequency span, and more. If None,
  523. no events are plotted.
  524. timeseries : array, shape (n_epochs, n_times) | None
  525. The stack of timeseries associated with the TFR epochs/trials.
  526. example_epochs : list | None
  527. The epoch/trial indices that will be used to plot example spectrogram
  528. trials alonside the average spectrogram.
  529. vlim : list, shape (2,) | None
  530. If not None, sets the colorbar lower and upper bounds for the
  531. spectrogram across all axes of the returned figure. If None, each
  532. spectrogram (including example epochs) sets its vlim independently.
  533. show_events : boolean
  534. Whether or not to overlay spectral event locations on the average
  535. spectrogram.
  537. Return
  538. ------
  539. fig : Matplotlib.figure.Figure instance
  540. The Figure instance containing the average spectrogram.
  541. '''
  543. if len(tfr.shape) != 3:
  544. raise ValueError(f'tfr should be a 3D array of shape (n_epochs, '
  545. f'n_freqs, n_times), got {tfr.shape}.')
  547. if example_epochs is not None:
  548. trial_idx_set = set(range(tfr.shape[0]))
  549. trial_idx_subset = set(example_epochs)
  550. if trial_idx_subset.intersection(trial_idx_set) != trial_idx_subset:
  551. raise ValueError('One or more of the specified example trial '
  552. 'indices does not exist in the provided tfr '
  553. 'array.')
  554. else:
  555. # set to empty list
  556. example_epochs = list()
  558. if show_events:
  559. spec_events_agg = sum(spec_events, [])
  560. else:
  561. spec_events_agg = None
  563. fig, axs = plt.subplots(nrows=len(example_epochs) + 1, ncols=1,
  564. sharex=True)
  566. # if example_epochs is None, ensure that axs is still subscriptable
  567. axs = np.atleast_1d(axs)
  569. # plot trial-average tfr
  570. tfr_avg = np.mean(tfr, axis=0).squeeze()
  571. plot_events(tfr=tfr_avg, times=times, freqs=freqs,
  572. event_band=event_band, spec_events=spec_events_agg,
  573. ax=axs[0], vlim=vlim, label='epoch avg.')
  575. # plot tfr + events for example trials
  576. if timeseries is not None and example_epochs is not None:
  577. max_ts_amplitude = np.max(timeseries[example_epochs])
  578. min_ts_amplitude = np.min(timeseries[example_epochs])
  579. ylim_ts = [max_ts_amplitude, min_ts_amplitude]
  581. for plot_idx, trial_idx in enumerate(example_epochs):
  582. # get spectral events for the current trial
  583. if spec_events is not None:
  584. trial_events = spec_events[trial_idx]
  585. else:
  586. trial_events = None
  588. # plot trial tfr
  589. tfr_trial = tfr[trial_idx, :, :].squeeze()
  590. timeseries_trial = timeseries[trial_idx, :]
  592. plot_events(tfr=tfr_trial, times=times, freqs=freqs,
  593. event_band=event_band, spec_events=trial_events,
  594. timeseries=timeseries_trial, ax=axs[plot_idx + 1],
  595. vlim=vlim, ylim_ts=ylim_ts, label=f'epoch {trial_idx}')
  597. axs[-1].set_xlabel('time (s)')
  598. axs[0].set_ylabel('freq. (Hz)')
  599. fig.tight_layout()
  601. return fig

spectralevents.py at commit cd0c83d, under BSD-3-Clause · at the source

Overview

Authors: Fumiaki Iwane1, William Hayward1, I M Dushyanthi Karunathilake1, Ethan R Buch1, Leonardo G Cohen1,2
ORCID iDs: Fumiaki Iwane, Leonardo G Cohen
  1. Human Cortical Physiology and Neurorehabilitation Section, National Institute of Neurological Disorders and Stroke, NIH, 9000 Rockville Pike, Bethesda, MD 20817, USA
  2. Lead contact
Institutions: National Institutes of Health (United States); National Institute of Neurological Disorders and Stroke (United States)
Journal: Cell reports, volume 45, issue 8, article 117845
Dates: published online 14 August 2026; in print August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.celrep.2026.117845 · PMID 42599800 · PMCID PMC13573638 · OpenAlex W7203483517
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: MEG (modality), cognitive (subfield)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging, Physiology & signal measures, Connectivity
Keywords: Memory, Hippocampus, Magnetoencephalography, MEG, Fatigue, Consolidation, Skill Learning, Motor Learning, Reactive Inhibition, Theta-gamma Phase-amplitude Coupling, Beta Oscillatory Activity, Cp: Neuroscience
Topic: Motor Control and Adaptation (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Intramural NIH HHS (Z01 NS003030, Z01 NS002978, ZIA NS002978, ZIA NS003030); National Institutes of Health
Citations: not cited yet (Europe PMC); 115 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

jonescompneurolab/SpectralEvents

License: BSD-3-Clause
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: cd0c83d2492446e03f1e58b79ed31bc73c024f22, 22 July 2024
Languages: MATLAB (6), Python (4), Jupyter (1)
Size: 69 files, 11 scripts
Software Heritage: not archived
Found in: the text, “MEG data analysis”
Holds: README, license file, environment (requirements.txt, setup.cfg), tests, continuous integration, 1 notebook
Not found: CITATION.cff, documentation
Tools: NumPy (3 files), SciPy (3 files), Matplotlib (2 files), Image Processing Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
13 files
At the source: github.com/jonescompneurolab/SpectralEvents

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;
  • 11 scripts, each with its path and the digest of its content;
  • 2 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

Datasets cited

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1016/j.celrep.2026.117845.

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 2, 28 September 2026

  • Publisher: n/a → Cell Press
  • Authors: added Leonardo G Cohen (0000-0002-1705-8773); removed Leonardo G Cohen

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 12 keywords, 2 funders, 113 references.

Cite

This paper

Iwane, F., Hayward, W., Karunathilake, I. M. D., Buch, E. R., & Cohen, L. G. (2026). Hippocampal skill memory expansion drives online performance dynamics during skill learning. Cell reports, 45(8), 117845. https://doi.org/10.1016/j.celrep.2026.117845

BibTeX

@article{iwane2026hippocampal,
author = {Iwane, Fumiaki and Hayward, William and Karunathilake, I M Dushyanthi and Buch, Ethan R and Cohen, Leonardo G},
title = {{Hippocampal skill memory expansion drives online performance dynamics during skill learning}},
journal = {Cell reports},
year = {2026},
month = aug,
volume = {45},
number = {8},
pages = {117845},
publisher = {Cell Press},
issn = {2211-1247},
doi = {10.1016/j.celrep.2026.117845},
url = {https://doi.org/10.1016/j.celrep.2026.117845},
pmid = {42599800},
pmcid = {PMC13573638}
}

RIS

TY - JOUR
AU - Iwane, Fumiaki
AU - Hayward, William
AU - Karunathilake, I M Dushyanthi
AU - Buch, Ethan R
AU - Cohen, Leonardo G
TI - Hippocampal skill memory expansion drives online performance dynamics during skill learning
T2 - Cell reports
J2 - Cell Rep
PY - 2026
DA - 2026/08/14
VL - 45
IS - 8
SP - 117845
SN - 2211-1247
PB - Cell Press
DO - 10.1016/j.celrep.2026.117845
UR - https://doi.org/10.1016/j.celrep.2026.117845
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.celrep.2026.117845",
"type": "article-journal",
"title": "Hippocampal skill memory expansion drives online performance dynamics during skill learning",
"container-title": "Cell reports",
"author": [
{
"family": "Iwane",
"given": "Fumiaki"
},
{
"family": "Hayward",
"given": "William"
},
{
"family": "Karunathilake",
"given": "I M Dushyanthi"
},
{
"family": "Buch",
"given": "Ethan R"
},
{
"family": "Cohen",
"given": "Leonardo G"
}
],
"container-title-short": "Cell Rep",
"volume": "45",
"issue": "8",
"page": "117845",
"DOI": "10.1016/j.celrep.2026.117845",
"PMID": "42599800",
"PMCID": "PMC13573638",
"ISSN": "2211-1247",
"publisher": "Cell Press",
"URL": "https://doi.org/10.1016/j.celrep.2026.117845",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
14
]
]
}
}

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/s41467-026-75345-6 [code]
Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval.
Journal: Nature communications
In common: Image Processing Toolbox, cognitive, 9 references
[2] doi:10.1016/j.isci.2026.116586 [code]
Condition-specific neural signatures of reactivation during post-retrieval rest: An EEG study.
Journal: iScience
In common: Image Processing Toolbox, SciPy, Matplotlib, 1 other tool, cognitive, 6 references
[3] doi:10.1038/s41593-026-02362-5 [code]
Replay of procedural memory is independent of the hippocampus.
Journal: Nature neuroscience
In common: Image Processing Toolbox, SciPy, Matplotlib, 1 other tool, cognitive, 4 references
[4] doi:10.1002/hbm.70562 [code]
Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning.
Journal: Human brain mapping
In common: SciPy, Matplotlib, NumPy, 5 references
[5] doi:10.1038/s41593-026-02285-1 [code]
Fixation duration on natural scenes is explained by memory encoding not processing demand.
Journal: Nature neuroscience
In common: SciPy, Matplotlib, NumPy, MEG, cognitive, 5 references
[6] doi:10.1162/imag.a.1203 [code]
Motor cortical areas facilitate schema-mediated integration of new motor information into memory.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: cognitive, 6 references
[7] doi:10.1038/s42003-026-10427-1 [code]
Phase-tuned modulation during reward expectancy in human anterior insular cortex.
Journal: Communications biology
In common: SciPy, NumPy, 6 references
[8] doi:10.1038/s41467-026-73818-2 [code]
Prefrontal parvalbumin neurons mediate working memory in a task demand-dependent manner.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, cognitive, 4 references
[9] doi:10.1162/imag.a.1359 [code]
Distinct roles of brain network flexibility in motor learning across age.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: 6 references
[10] doi:10.1002/hbm.70579
Neural Correlates of Motor Sequence Learning and Enhanced Offline Consolidation in 7-11-Year-Old Children.
Journal: Human brain mapping
In common: cognitive, 5 references

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.