OSCR

Transient gamma events delineate somatosensory modality in S1.

Code ↔ Paper

3 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 3 matches
  1. [1] § Results › Transient gamma events correlate with stimulus modality and nocifensive response ↔ spectralevents.py, lines 293–399 · score 0.76 · un normalized TFR, event frequency span, Event onset, suprathreshold regions, maximum power, spectral events
  2. [2] § Results › Transient gamma events correlate with stimulus modality and nocifensive response ↔ spectralevents_find.m, lines 282–404 · score 0.71 · un normalized TFR, Event onset, suprathreshold regions, maximum power, spectral events, magnitude
  3. [3] § STAR★Methods › Method details › Spectral event analysis ↔ spectralevents_find.m, lines 203–280 · score 0.68 · spectralEvents, event metrics, un normalized, Shin, frequency bands, peaks

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

MATLAB · 528 lines · 32 KB · BSD-3-Clause · 2 matches

  1. function specEv_struct = spectralevents_find(findMethod, eventBand, thrFOM, tVec, fVec, TFR, classLabels)
  2. % SPECTRALEVENTS_FIND Algorithm for finding and calculating spectral
  3. % events on a trial-by-trial basis of a single subject/session. Uses
  4. % one of three methods before further analyzing and organizing event
  5. % features:
  6. %
  7. % 1) (Primary event detection method in Shin et al. eLife 2017): Find
  8. % spectral events by first retrieving all local maxima in
  9. % un-normalized TFR using imregionalmax, then selecting suprathreshold
  10. % peaks within the frequency band of interest. This method allows for
  11. % multiple, overlapping events to occur in a given suprathreshold
  12. % region and does not guarantee the presence of within-band,
  13. % suprathreshold activity in any given trial will render an event.
  14. % 2) Find spectral events by first thresholding
  15. % entire normalize TFR (over all frequencies), then finding local
  16. % maxima. Discard those of lesser magnitude in each suprathreshold
  17. % region, respectively, s.t. only the greatest local maximum in each
  18. % region survives (when more than one local maxima in a region have
  19. % the same greatest value, their respective event timing, freq.
  20. % location, and boundaries at full-width half-max are calculated
  21. % separately and averaged). This method does not allow for overlapping
  22. % events to occur in a given suprathreshold region and does not
  23. % guarantee the presence of within-band, suprathreshold activity in
  24. % any given trial will render an event.
  25. % 3) Find spectral events by first thresholding
  26. % normalized TFR in frequency band of interest, then finding local
  27. % maxima. Discard those of lesser magnitude in each suprathreshold region,
  28. % respectively, s.t. only the greatest local maximum in each region
  29. % survives (when more than one local maxima in a region have the same
  30. % greatest value, their respective event timing, freq. location, and
  31. % boundaries at full-width half-max are calculated separately and
  32. % averaged). This method does not allow for overlapping events to occur in
  33. % a given suprathreshold region and ensures the presence of
  34. % within-band, suprathreshold activity in any given trial will render
  35. % an event.
  36. %
  37. % specEv_struct = SPECTRALEVENTS_FIND(findMethod,eventBand,thrFOM,tVec,fVec,TFR,classLabels)
  38. %
  39. % Inputs:
  40. % findMethod - integer value specifying which event-finding method to use
  41. % (1, 2, or 3). Note that the method specifies how much overlap
  42. % exists between events. Use 1 to replicate the method used in
  43. % et al. eLife 2017.
  44. % eventBand - range of frequencies ([Fmin_event Fmax_event]; Hz) over
  45. % which above-threshold spectral power events are classified.
  46. % thrFOM - factors of median threshold; positive real number used to
  47. % threshold local maxima and classify events (see Shin et al. eLife
  48. % 2017 for discussion concerning this value).
  49. % tVec - time vector (s) over which the time-frequency response (TFR) is
  50. % calcuated.
  51. % fVec - frequency vector (Hz) over which the time-frequency response
  52. % (TFR) is calcuated.
  53. % TFR - time-frequency response (TFR) (frequency-by-time-trial) for a
  54. % single subject/session.
  55. % classLabels - numeric or logical 1-row array of trial classification
  56. % labels; associates each trial of the given subject/session to an
  57. % experimental condition/outcome/state (e.g., hit or miss, detect or
  58. % non-detect, attend-to or attend away).
  59. %
  60. % Outputs:
  61. % specEv_struct - event feature structure with three main sub-structures:
  62. % TrialSummary (trial-level features), Events (individual event
  63. % characteristics), and IEI (inter-event intervals from all trials
  64. % and those associated with only a given class label).
  65. %
  66. % See also SPECTRALEVENTS, SPECTRALEVENTS_FIND, SPECTRALEVENTS_TS2TFR, SPECTRALEVENTS_VIS.
  68. % Initialize general data parameters
  69. eventBand_inds = fVec>=eventBand(1) & fVec<=eventBand(2); %Logical vector representing indices of freq vector within eventBand
  70. if size(eventBand_inds,1)~=length(eventBand_inds)
  71. eventBand_inds = eventBand_inds'; %Transpose so that the dimensions correspond with the frequency-domain dimension of the TFR
  72. end
  73. flength = size(TFR,1); %Number of elements in discrete frequency spectrum
  74. tlength = size(TFR,2); %Number of points in time
  75. numTrials = size(TFR,3); %Number of trials
  76. classes = unique(classLabels);
  77. medianpower = median(reshape(TFR, size(TFR,1), size(TFR,2)*size(TFR,3)), 2); %Median power at each frequency across all trials
  78. thr = thrFOM*medianpower; %Spectral event threshold for each frequency value
  80. % Validate consistency of parameter dimensions
  81. if flength~=length(fVec) || tlength~=length(tVec) || numTrials~=length(classLabels)
  82. error('Mismatch in input parameter dimensions!')
  83. end
  85. % Find events using the method-of-choice
  86. spectralEvents = []; %Array for storing event results
  87. switch findMethod
  88. case 1
  89. find_localmax_method_1;
  90. case 2
  91. find_localmax_method_2;
  92. case 3
  93. find_localmax_method_3;
  94. otherwise
  95. error('Unknown event-finding method.')
  96. end
  98. % Make sure this subject/session contains >1 events
  99. if isempty(spectralEvents)
  100. disp('Warning!! This subject/session contains no events!!')
  101. specEv_struct = struct('TrialSummary',[],'Events',[],'IEI',[]);
  102. return;
  103. end
  105. % Identify and organize event features
  107. % Matrix of event features: each row is an event
  108. % 11 column matrix with 1. trial index, 2. hit/miss, 3. maxima frequency, 4. lowerbound frequency, 5. upperbound frequency, 6. frequency span, ...
  109. % 7. maxima timing, 8. event onset timing, 9. event offset timing, 10. event duration, 11. maxima power, 12. maxima/median power
  110. spectralEvents_columnlabel={'trialind', 'classLabels', 'maximafreq', 'lowerboundFspan', 'upperboundFspan', 'Fspan', ...
  111. 'maximatiming', 'onsettiming', 'offsettiming', 'duration', 'maximapower', 'maximapowerFOM'};
  112. for rci=1:numel(spectralEvents_columnlabel)
  113. eventsind.(spectralEvents_columnlabel{rci})=rci;
  114. end
  116. trialSummary.classLabels = classLabels';
  118. trialSummary.meanpower = mean(squeeze(mean(TFR(eventBand_inds,:,:),2)) ./ repmat(medianpower(eventBand_inds),1,numTrials), 1)'; %Mean trial power normalized to frequency-specific median
  120. suprathrTFR = TFR>=repmat(thr,1,tlength,numTrials);
  121. trialSummary.coverage = squeeze(sum(sum(suprathrTFR(eventBand_inds,:,:),1),2)) *100 / (nnz(eventBand_inds)*tlength); %Calculated in percentage
  123. % Initialize column vectors
  124. trialSummary.eventnumber = nan(numTrials,1);
  125. trialSummary.meaneventpower = nan(numTrials,1);
  126. trialSummary.meaneventduration = nan(numTrials,1);
  127. trialSummary.meaneventFspan = nan(numTrials,1);
  128. trialSummary.mostrecenteventtiming = nan(numTrials,1);
  129. trialSummary.mostrecenteventpower = nan(numTrials,1);
  130. trialSummary.mostrecenteventduration = nan(numTrials,1);
  131. trialSummary.mostrecenteventFspan = nan(numTrials,1);
  133. % Iterate through trials
  134. for tri=1:numTrials
  135. trialSummary.eventnumber(tri)=nnz(spectralEvents(:,1)==tri);
  136. if nnz(spectralEvents(:,1)==tri)==0
  137. trialSummary.meaneventpower(tri) = 0; % traces2TFR always returns a positive value
  138. trialSummary.meaneventduration(tri) = 0;
  139. trialSummary.meaneventFspan(tri) = 0;
  140. trialSummary.mostrecenteventtiming(tri) = tVec(1)-mean(diff(tVec));
  141. trialSummary.mostrecenteventpower(tri) = 0;
  142. trialSummary.mostrecenteventduration(tri) = 0;
  143. trialSummary.mostrecenteventFspan(tri) = 0;
  144. else
  145. trialSummary.meaneventpower(tri) = mean(spectralEvents(spectralEvents(:,eventsind.trialind)==tri,eventsind.maximapowerFOM)); % traces2TFR always returns a positive value
  146. trialSummary.meaneventduration(tri) = mean(spectralEvents(spectralEvents(:,eventsind.trialind)==tri,eventsind.duration));
  147. trialSummary.meaneventFspan(tri) = mean(spectralEvents(spectralEvents(:,eventsind.trialind)==tri,eventsind.Fspan));
  148. trialSummary.mostrecenteventtiming(tri) = spectralEvents(find(spectralEvents(:,eventsind.trialind)==tri,1,'last'), eventsind.maximatiming);
  149. trialSummary.mostrecenteventpower(tri) = spectralEvents(find(spectralEvents(:,eventsind.trialind)==tri,1,'last'), eventsind.maximapowerFOM);
  150. trialSummary.mostrecenteventduration(tri) = spectralEvents(find(spectralEvents(:,eventsind.trialind)==tri,1,'last'), eventsind.duration);
  151. trialSummary.mostrecenteventFspan(tri) = spectralEvents(find(spectralEvents(:,eventsind.trialind)==tri,1,'last'), eventsind.Fspan);
  152. end
  153. end
  155. % Event dependent features (mean power, mean length, most recent timing):
  156. % need special treatment for zero event trials
  157. specialFeat.field = {'meaneventpower','meaneventduration','meaneventFspan','mostrecenteventtiming',...
  158. 'mostrecenteventpower','mostrecenteventduration','mostrecenteventFspan'};
  160. % Percent change from mean (PCM)
  161. trialSum_featNames = fieldnames(trialSummary);
  162. for feat_i=2:numel(trialSum_featNames)
  163. pcm_name = [trialSum_featNames{feat_i},'_pcm'];
  164. feature = trialSummary.(trialSum_featNames{feat_i});
  166. % Control for features that need special treatment
  167. validtrials = trialSummary.eventnumber>0 ; %Trials that do have events
  168. if ismember(trialSum_featNames{feat_i},specialFeat.field)
  169. trialSummary.(pcm_name) = 100 * (feature-mean(feature(validtrials))) ./ repmat(abs(mean(feature(validtrials))),numTrials,1);
  170. else
  171. trialSummary.(pcm_name) = 100 * (feature-mean(feature))./repmat(abs(mean(feature)),numTrials,1);
  172. end
  173. end
  175. % Inter-event interval (IEI)
  176. ieitemp=diff(spectralEvents(:,eventsind.maximatiming));
  177. sametrial=(diff(spectralEvents(:,eventsind.trialind))==0);
  178. IEI.IEI_all = ieitemp(sametrial);
  180. for cls_i=1:numel(classes)
  181. fieldName = ['IEI_',num2str(classes(cls_i))];
  182. iei_class=diff(spectralEvents(spectralEvents(:,eventsind.classLabels)==classes(cls_i),eventsind.maximatiming));
  183. sametrial_class=(diff(spectralEvents(spectralEvents(:,eventsind.classLabels)==classes(cls_i),eventsind.trialind)) == 0);
  184. IEI.(fieldName) = iei_class(sametrial_class);
  185. end
  187. % Assign output structure with 3 main branches: trial-level summary
  188. % (TrialSummary), trial-specific events (Events), and mean inter-event
  189. % interval across trials (IEI)
  190. specEv_struct.TrialSummary = struct('NumTrials',numTrials,'SpecialFeatures',specialFeat,'TrialSummary',trialSummary);
  191. specEv_struct.Events = struct('EventBand',eventBand,'ThrFOM',thrFOM,'MedianPower',medianpower,'Threshold',thr,'Events',struct(spectralEvents_columnlabel{1},spectralEvents(:,1),...
  192. spectralEvents_columnlabel{2},spectralEvents(:,2),spectralEvents_columnlabel{3},spectralEvents(:,3),spectralEvents_columnlabel{4},spectralEvents(:,4),...
  193. spectralEvents_columnlabel{5},spectralEvents(:,5),spectralEvents_columnlabel{6},spectralEvents(:,6),spectralEvents_columnlabel{7},spectralEvents(:,7),...
  194. spectralEvents_columnlabel{8},spectralEvents(:,8),spectralEvents_columnlabel{9},spectralEvents(:,9),spectralEvents_columnlabel{10},spectralEvents(:,10),...
  195. spectralEvents_columnlabel{11},spectralEvents(:,11),spectralEvents_columnlabel{12},spectralEvents(:,12)));
  196. specEv_struct.IEI = IEI;
  203. function find_localmax_method_1
  204. % 1st event-finding method (primary event detection method in Shin et
  205. % al. eLife 2017): Find spectral events by first retrieving all local
  206. % maxima in un-normalized TFR using imregionalmax, then selecting
  207. % suprathreshold peaks within the frequency band of interest. This
  208. % method allows for multiple, overlapping events to occur in a given
  209. % suprathreshold region and does not guarantee the presence of
  210. % within-band, suprathreshold activity in any given trial will render
  211. % an event.
  213. % spectralEvents: 12 column matrix for storing local max event metrics: trial
  214. % index, hit/miss, maxima frequency, lowerbound frequency, upperbound
  215. % frequency, frequency span, maxima timing, event onset timing, event
  216. % offset timing, event duration, maxima power, maxima/median power
  217. spectralEvents = [];
  219. % Finds_localmax: stores peak frequecy at each local max (columns) for each
  220. % trial (rows)
  221. Finds_localmax = [];
  223. % Retrieve all local maxima in TFR using imregionalmax
  224. for ti=1:numTrials
  225. [peakF,peakT] = find(imregionalmax(squeeze(TFR(:,:,ti)))); %Indices of max local power
  226. peakpower = TFR(find(imregionalmax(squeeze(TFR(:,:,ti))))+(ti-1)*flength*tlength); %Power values at local maxima (vector; compiles across frequencies and time)
  228. % Find local maxima lowerbound, upperbound, and full width at half max
  229. % for both frequency and time
  230. Ffwhm = NaN(numel(peakpower),3); %2D matrix for freq-dimension event metrics with columns containing lowerbound, upperbound, and fwhm, respectively
  231. Tfwhm = NaN(numel(peakpower),3); %2D matrix for time-dimension event metrics with columns containing lowerbound, upperbound, and fwhm, respectively
  232. for lmi=1:numel(peakpower)
  233. lmF_underthr = find(squeeze(TFR(:,peakT(lmi),ti) < peakpower(lmi)/2)); %Indices of TFR frequencies of < half max power at the time of a given local peak
  234. if ~isempty(find(lmF_underthr < peakF(lmi), 1)) && ~isempty(find(lmF_underthr > peakF(lmi), 1))
  235. Ffwhm(lmi,1) = fVec(lmF_underthr(find(lmF_underthr < peakF(lmi),1,'last'))+1);
  236. Ffwhm(lmi,2) = fVec(lmF_underthr(find(lmF_underthr > peakF(lmi),1,'first'))-1);
  237. Ffwhm(lmi,3) = Ffwhm(lmi,2)-Ffwhm(lmi,1)+ min(diff(fVec));
  238. elseif isempty(find(lmF_underthr < peakF(lmi),1)) && ~isempty(find(lmF_underthr > peakF(lmi),1))
  239. Ffwhm(lmi,1) = fVec(1);
  240. Ffwhm(lmi,2) = fVec(lmF_underthr(find(lmF_underthr > peakF(lmi),1,'first'))-1);
  241. Ffwhm(lmi,3) = 2*(Ffwhm(lmi,2)-fVec(peakF(lmi)))+ min(diff(fVec));
  242. elseif ~isempty(find(lmF_underthr < peakF(lmi),1)) && isempty(find(lmF_underthr > peakF(lmi),1))
  243. Ffwhm(lmi,1) = fVec(lmF_underthr(find(lmF_underthr < peakF(lmi),1,'last'))+1);
  244. Ffwhm(lmi,2) = fVec(end);
  245. Ffwhm(lmi,3) = 2*(fVec(peakF(lmi))-Ffwhm(lmi,1))+ min(diff(fVec));
  246. else
  247. Ffwhm(lmi,1) = fVec(1);
  248. Ffwhm(lmi,2) = fVec(end);
  249. Ffwhm(lmi,3) = 2*(fVec(end)-fVec(1)+min(diff(fVec)));
  250. end
  252. lmT_underthr = find(squeeze(TFR(peakF(lmi),:,ti) < peakpower(lmi)/2)); %Indices of TFR times of < half max power at the freq of a given local peak
  253. if ~isempty(find(lmT_underthr < peakT(lmi), 1)) && ~isempty(find(lmT_underthr > peakT(lmi), 1))
  254. Tfwhm(lmi,1) = tVec(lmT_underthr(find(lmT_underthr < peakT(lmi),1,'last'))+1);
  255. Tfwhm(lmi,2) = tVec(lmT_underthr(find(lmT_underthr > peakT(lmi),1,'first'))-1);
  256. Tfwhm(lmi,3) = Tfwhm(lmi,2)-Tfwhm(lmi,1)+ min(diff(tVec));
  257. elseif isempty(find(lmT_underthr < peakT(lmi),1)) && ~isempty(find(lmT_underthr > peakT(lmi),1))
  258. Tfwhm(lmi,1) = tVec(1);
  259. Tfwhm(lmi,2) = tVec(lmT_underthr(find(lmT_underthr > peakT(lmi),1,'first'))-1);
  260. Tfwhm(lmi,3) = 2*(Tfwhm(lmi,2)-tVec(peakT(lmi)))+ min(diff(tVec));
  261. elseif ~isempty(find(lmT_underthr < peakT(lmi),1)) && isempty(find(lmT_underthr > peakT(lmi),1))
  262. Tfwhm(lmi,1) = tVec(lmT_underthr(find(lmT_underthr < peakT(lmi),1,'last'))+1);
  263. Tfwhm(lmi,2) = tVec(end);
  264. Tfwhm(lmi,3) = 2*(tVec(peakT(lmi))-Tfwhm(lmi,1))+ min(diff(tVec));
  265. else
  266. Tfwhm(lmi,1) = tVec(1);
  267. Tfwhm(lmi,2) = tVec(end);
  268. Tfwhm(lmi,3) = 2*(tVec(end)-tVec(1)+min(diff(tVec)));
  269. end
  270. end
  272. % 12 column matrix with 1. trial index, 2. trial class, 3. maxima frequency, 4. lowerbound frequency, 5. upperbound frequency, 6. frequency span, ...
  273. % 7. maxima timing, 8. event onset timing, 9. event offset timing, 10. event duration, 11. maxima power, 12. maxima/median power, ...
  274. spectralEvents = [spectralEvents; ti*ones(size(peakF)) classLabels(ti)*ones(size(peakF)) fVec(peakF)' Ffwhm tVec(peakT)' Tfwhm peakpower peakpower./medianpower(peakF)];
  275. Finds_localmax = [Finds_localmax; peakF];
  276. end
  278. % Pick out maxima above threshold and within the frequency band of interest
  279. spectralEvents = spectralEvents((spectralEvents(:,3)>=eventBand(1) & spectralEvents(:,3)<=eventBand(2) & spectralEvents(:,11)>=thr(Finds_localmax)),:); %Select local maxima
  280. end
  282. function find_localmax_method_2
  283. % 2nd event-finding method: Find spectral events by first thresholding
  284. % entire normalize TFR (over all frequencies), then finding local
  285. % maxima. This method does not allow for overlapping events to occur in
  286. % a given suprathreshold region and does not guarantee the presence of
  287. % within-band, suprathreshold activity in any given trial will render
  288. % an event.
  290. % spectralEvents: 12 column matrix for storing local max event metrics: trial
  291. % index, hit/miss, maxima frequency, lowerbound frequency, upperbound
  292. % frequency, frequency span, maxima timing, event onset timing, event
  293. % offset timing, event duration, maxima power, maxima/median power
  294. spectralEvents = [];
  296. % Retrieve local maxima in normalized TFR using imregionalmax,
  297. % discard those of lesser (un-normalized) magnitude in each suprathreshold
  298. % region, respectively, and characterize event boundaries (at half max)
  299. for ti=1:numTrials
  300. TFR_ST = squeeze(TFR(:,:,ti))./medianpower; %Suprathreshold TFR: first isolate 2D TFR matrix and normalize
  301. TFR_ST(TFR_ST<thrFOM) = 0; %Set infrathreshold values to zero
  303. % Find all local maxima in suprathreshold TFR
  304. TFR_LM = TFR_ST.*imregionalmax(TFR_ST); %Threshold TFR at each respective local maximum
  305. numTotalPeaks = nnz(TFR_LM);
  307. % Escape this iteration when this trial contains no suprathreshold
  308. % local maxima
  309. if numTotalPeaks==0
  310. continue
  311. end
  313. % Find max peak in each respective suprathreshold region
  314. [~,regions,numReg,~] = bwboundaries(TFR_ST>=thrFOM); %Separate suprathreshold regions
  315. evPeakF = cell(1,numReg);
  316. evPeakT = cell(1,numReg);
  317. evPeakpower = nan(numReg,1);
  318. for reg_i=1:numReg
  319. region = zeros(size(TFR_ST)); %Initialize a blank image that will contain a single region
  320. region(regions==reg_i) = 1; %Set elements (pixels) in region to the value 1
  321. TFR_reg = TFR_LM.*region; %Regional local maxima
  322. [peakF_reg,peakT_reg] = find(TFR_reg); %Indices of regional local maxima
  323. peakpower_reg = TFR(find(TFR_reg)+(ti-1)*flength*tlength); %Power values at regional local maxima
  324. maxPeakpower = max(peakpower_reg);
  325. maxPeak_inds = find(peakpower_reg==maxPeakpower); %Indices of all instances where local maxima have the max peak power
  326. evPeakF{reg_i} = peakF_reg(maxPeak_inds); %Select TFR indices at max regional peak
  327. evPeakT{reg_i} = peakT_reg(maxPeak_inds); %Select TFR indices at max regional peak
  328. evPeakpower(reg_i) = maxPeakpower(1);
  329. end
  331. % Find local maxima lowerbound, upperbound, and full width at half max
  332. % for both frequency and time
  333. evBndsF = nan(numReg,3);
  334. evBndsT = nan(numReg,3);
  335. evPeakF_inds = nan(numReg,1);
  336. evPeakT_inds = nan(numReg,1);
  337. evPeakpower_norm = nan(numReg,1);
  338. for reg_i=1:numReg
  339. numRegPeaks = numel(evPeakF{reg_i});
  340. peakF = evPeakF{reg_i};
  341. peakT = evPeakT{reg_i};
  342. peakpower = evPeakpower(reg_i);
  344. Ffwhm = nan(numRegPeaks,3); %2D matrix for freq-dimension event metrics with columns containing lowerbound, upperbound, and fwhm, respectively
  345. Tfwhm = nan(numRegPeaks,3); %2D matrix for time-dimension event metrics with columns containing lowerbound, upperbound, and fwhm, respectively
  346. peakpower_norm = nan(numRegPeaks,1); %Vector for storing the normalized power at each regional peak
  347. for lmi=1:numRegPeaks
  348. lmF_underthr = find(squeeze(TFR(:,peakT(lmi),ti) < peakpower/2)); %Indices of TFR frequencies of < half max power at the time of a given local peak
  349. if ~isempty(find(lmF_underthr < peakF(lmi), 1)) && ~isempty(find(lmF_underthr > peakF(lmi), 1))
  350. Ffwhm(lmi,1) = fVec(lmF_underthr(find(lmF_underthr < peakF(lmi),1,'last'))+1);
  351. Ffwhm(lmi,2) = fVec(lmF_underthr(find(lmF_underthr > peakF(lmi),1,'first'))-1);
  352. Ffwhm(lmi,3) = Ffwhm(lmi,2)-Ffwhm(lmi,1)+ min(diff(fVec));
  353. elseif isempty(find(lmF_underthr < peakF(lmi),1)) && ~isempty(find(lmF_underthr > peakF(lmi),1))
  354. Ffwhm(lmi,1) = fVec(1);
  355. Ffwhm(lmi,2) = fVec(lmF_underthr(find(lmF_underthr > peakF(lmi),1,'first'))-1);
  356. Ffwhm(lmi,3) = 2*(Ffwhm(lmi,2)-fVec(peakF(lmi)))+ min(diff(fVec));
  357. elseif ~isempty(find(lmF_underthr < peakF(lmi),1)) && isempty(find(lmF_underthr > peakF(lmi),1))
  358. Ffwhm(lmi,1) = fVec(lmF_underthr(find(lmF_underthr < peakF(lmi),1,'last'))+1);
  359. Ffwhm(lmi,2) = fVec(end);
  360. Ffwhm(lmi,3) = 2*(fVec(peakF(lmi))-Ffwhm(lmi,1))+ min(diff(fVec));
  361. else
  362. Ffwhm(lmi,1) = fVec(1);
  363. Ffwhm(lmi,2) = fVec(end);
  364. Ffwhm(lmi,3) = 2*(fVec(end)-fVec(1)+min(diff(fVec)));
  365. end
  367. lmT_underthr = find(squeeze(TFR(peakF(lmi),:,ti) < peakpower/2)); %Indices of TFR times of < half max power at the freq of a given local peak
  368. if ~isempty(find(lmT_underthr < peakT(lmi), 1)) && ~isempty(find(lmT_underthr > peakT(lmi), 1))
  369. Tfwhm(lmi,1) = tVec(lmT_underthr(find(lmT_underthr < peakT(lmi),1,'last'))+1);
  370. Tfwhm(lmi,2) = tVec(lmT_underthr(find(lmT_underthr > peakT(lmi),1,'first'))-1);
  371. Tfwhm(lmi,3) = Tfwhm(lmi,2)-Tfwhm(lmi,1)+ min(diff(tVec));
  372. elseif isempty(find(lmT_underthr < peakT(lmi),1)) && ~isempty(find(lmT_underthr > peakT(lmi),1))
  373. Tfwhm(lmi,1) = tVec(1);
  374. Tfwhm(lmi,2) = tVec(lmT_underthr(find(lmT_underthr > peakT(lmi),1,'first'))-1);
  375. Tfwhm(lmi,3) = 2*(Tfwhm(lmi,2)-tVec(peakT(lmi)))+ min(diff(tVec));
  376. elseif ~isempty(find(lmT_underthr < peakT(lmi),1)) && isempty(find(lmT_underthr > peakT(lmi),1))
  377. Tfwhm(lmi,1) = tVec(lmT_underthr(find(lmT_underthr < peakT(lmi),1,'last'))+1);
  378. Tfwhm(lmi,2) = tVec(end);
  379. Tfwhm(lmi,3) = 2*(tVec(peakT(lmi))-Tfwhm(lmi,1))+ min(diff(tVec));
  380. else
  381. Tfwhm(lmi,1) = tVec(1);
  382. Tfwhm(lmi,2) = tVec(end);
  383. Tfwhm(lmi,3) = 2*(tVec(end)-tVec(1)+min(diff(tVec)));
  384. end
  386. peakpower_norm(lmi) = TFR_ST(peakF(lmi),peakT(lmi));
  387. end
  389. evBndsF(reg_i,:) = mean(Ffwhm,1);
  390. evBndsT(reg_i,:) = mean(Tfwhm,1);
  391. evPeakF_inds(reg_i) = round(mean(peakF));
  392. evPeakT_inds(reg_i) = round(mean(peakT));
  393. evPeakpower_norm(reg_i) = mean(peakpower_norm);
  394. end
  396. % 12 column matrix with 1. trial index, 2. trial class, 3. maxima frequency, 4. lowerbound frequency, 5. upperbound frequency, 6. frequency span, ...
  397. % 7. maxima timing, 8. event onset timing, 9. event offset timing, 10. event duration, 11. maxima power, 12. maxima/median power, ...
  398. spectralEvents = [spectralEvents; ti*ones(size(evPeakpower)) classLabels(ti)*ones(size(evPeakpower))...
  399. fVec(evPeakF_inds)' evBndsF tVec(evPeakT_inds)' evBndsT evPeakpower evPeakpower_norm];
  400. end
  402. % Pick out maxima within the frequency band of interest
  403. spectralEvents = spectralEvents((spectralEvents(:,3)>=eventBand(1) & spectralEvents(:,3)<=eventBand(2)),:); %Select local maxima
  404. end
  406. function find_localmax_method_3
  407. % 3rd event-finding method: Find spectral events by first thresholding
  408. % normalized TFR in frequency band of interest, then finding local
  409. % maxima. This method does not allow for overlapping events to occur in
  410. % a given suprathreshold region and ensures the presence of
  411. % within-band, suprathreshold activity in any given trial will render
  412. % an event.
  414. % spectralEvents: 12 column matrix for storing local max event metrics: trial
  415. % index, hit/miss, maxima frequency, lowerbound frequency, upperbound
  416. % frequency, frequency span, maxima timing, event onset timing, event
  417. % offset timing, event duration, maxima power, maxima/median power
  418. spectralEvents = [];
  420. % Retrieve local maxima in normalized TFR using imregionalmax,
  421. % discard those of lesser (un-normalized) magnitude in each suprathreshold
  422. % region, respectively, and characterize event boundaries (at half max)
  423. for ti=1:numTrials
  424. TFR_ST = squeeze(TFR(:,:,ti))./medianpower; %Suprathreshold TFR: first isolate 2D TFR matrix and normalize
  425. TFR_ST(TFR_ST<thrFOM) = 0; %Set infrathreshold values to zero
  426. TFR_ST = TFR_ST.*eventBand_inds; %Set out-of-band values to zero
  428. % Find all local maxima in suprathreshold TFR
  429. TFR_LM = TFR_ST.*imregionalmax(TFR_ST); %Threshold TFR at each respective local maximum
  430. numTotalPeaks = nnz(TFR_LM);
  432. % Escape this iteration when this trial contains no suprathreshold
  433. % local maxima
  434. if numTotalPeaks==0
  435. continue
  436. end
  438. % Find max peak in each respective suprathreshold region
  439. [~,regions,numReg,~] = bwboundaries(TFR_ST>=thrFOM); %Separate suprathreshold regions
  440. evPeakF = cell(1,numReg);
  441. evPeakT = cell(1,numReg);
  442. evPeakpower = nan(numReg,1);
  443. for reg_i=1:numReg
  444. region = zeros(size(TFR_ST)); %Initialize a blank image that will contain a single region
  445. region(regions==reg_i) = 1; %Set elements (pixels) in region to the value 1
  446. TFR_reg = TFR_LM.*region; %Regional local maxima
  447. [peakF_reg,peakT_reg] = find(TFR_reg); %Indices of regional local maxima
  448. peakpower_reg = TFR(find(TFR_reg)+(ti-1)*flength*tlength); %Power values at regional local maxima
  449. maxPeakpower = max(peakpower_reg);
  450. maxPeak_inds = find(peakpower_reg==maxPeakpower); %Indices of all instances where local maxima have the max peak power
  451. evPeakF{reg_i} = peakF_reg(maxPeak_inds); %Select TFR indices at max regional peak
  452. evPeakT{reg_i} = peakT_reg(maxPeak_inds); %Select TFR indices at max regional peak
  453. evPeakpower(reg_i) = maxPeakpower(1);
  454. end
  456. % Find local maxima lowerbound, upperbound, and full width at half max
  457. % for both frequency and time
  458. evBndsF = nan(numReg,3);
  459. evBndsT = nan(numReg,3);
  460. evPeakF_inds = nan(numReg,1);
  461. evPeakT_inds = nan(numReg,1);
  462. evPeakpower_norm = nan(numReg,1);
  463. for reg_i=1:numReg
  464. numRegPeaks = numel(evPeakF{reg_i});
  465. peakF = evPeakF{reg_i};
  466. peakT = evPeakT{reg_i};
  467. peakpower = evPeakpower(reg_i);
  469. Ffwhm = nan(numRegPeaks,3); %2D matrix for freq-dimension event metrics with columns containing lowerbound, upperbound, and fwhm, respectively
  470. Tfwhm = nan(numRegPeaks,3); %2D matrix for time-dimension event metrics with columns containing lowerbound, upperbound, and fwhm, respectively
  471. peakpower_norm = nan(numRegPeaks,1); %Vector for storing the normalized power at each regional peak
  472. for lmi=1:numRegPeaks
  473. lmF_underthr = find(squeeze(TFR(:,peakT(lmi),ti) < peakpower/2)); %Indices of TFR frequencies of < half max power at the time of a given local peak
  474. if ~isempty(find(lmF_underthr < peakF(lmi), 1)) && ~isempty(find(lmF_underthr > peakF(lmi), 1))
  475. Ffwhm(lmi,1) = fVec(lmF_underthr(find(lmF_underthr < peakF(lmi),1,'last'))+1);
  476. Ffwhm(lmi,2) = fVec(lmF_underthr(find(lmF_underthr > peakF(lmi),1,'first'))-1);
  477. Ffwhm(lmi,3) = Ffwhm(lmi,2)-Ffwhm(lmi,1)+ min(diff(fVec));
  478. elseif isempty(find(lmF_underthr < peakF(lmi),1)) && ~isempty(find(lmF_underthr > peakF(lmi),1))
  479. Ffwhm(lmi,1) = fVec(1);
  480. Ffwhm(lmi,2) = fVec(lmF_underthr(find(lmF_underthr > peakF(lmi),1,'first'))-1);
  481. Ffwhm(lmi,3) = 2*(Ffwhm(lmi,2)-fVec(peakF(lmi)))+ min(diff(fVec));
  482. elseif ~isempty(find(lmF_underthr < peakF(lmi),1)) && isempty(find(lmF_underthr > peakF(lmi),1))
  483. Ffwhm(lmi,1) = fVec(lmF_underthr(find(lmF_underthr < peakF(lmi),1,'last'))+1);
  484. Ffwhm(lmi,2) = fVec(end);
  485. Ffwhm(lmi,3) = 2*(fVec(peakF(lmi))-Ffwhm(lmi,1))+ min(diff(fVec));
  486. else
  487. Ffwhm(lmi,1) = fVec(1);
  488. Ffwhm(lmi,2) = fVec(end);
  489. Ffwhm(lmi,3) = 2*(fVec(end)-fVec(1)+min(diff(fVec)));
  490. end
  492. lmT_underthr = find(squeeze(TFR(peakF(lmi),:,ti) < peakpower/2)); %Indices of TFR times of < half max power at the freq of a given local peak
  493. if ~isempty(find(lmT_underthr < peakT(lmi), 1)) && ~isempty(find(lmT_underthr > peakT(lmi), 1))
  494. Tfwhm(lmi,1) = tVec(lmT_underthr(find(lmT_underthr < peakT(lmi),1,'last'))+1);
  495. Tfwhm(lmi,2) = tVec(lmT_underthr(find(lmT_underthr > peakT(lmi),1,'first'))-1);
  496. Tfwhm(lmi,3) = Tfwhm(lmi,2)-Tfwhm(lmi,1)+ min(diff(tVec));
  497. elseif isempty(find(lmT_underthr < peakT(lmi),1)) && ~isempty(find(lmT_underthr > peakT(lmi),1))
  498. Tfwhm(lmi,1) = tVec(1);
  499. Tfwhm(lmi,2) = tVec(lmT_underthr(find(lmT_underthr > peakT(lmi),1,'first'))-1);
  500. Tfwhm(lmi,3) = 2*(Tfwhm(lmi,2)-tVec(peakT(lmi)))+ min(diff(tVec));
  501. elseif ~isempty(find(lmT_underthr < peakT(lmi),1)) && isempty(find(lmT_underthr > peakT(lmi),1))
  502. Tfwhm(lmi,1) = tVec(lmT_underthr(find(lmT_underthr < peakT(lmi),1,'last'))+1);
  503. Tfwhm(lmi,2) = tVec(end);
  504. Tfwhm(lmi,3) = 2*(tVec(peakT(lmi))-Tfwhm(lmi,1))+ min(diff(tVec));
  505. else
  506. Tfwhm(lmi,1) = tVec(1);
  507. Tfwhm(lmi,2) = tVec(end);
  508. Tfwhm(lmi,3) = 2*(tVec(end)-tVec(1)+min(diff(tVec)));
  509. end
  511. peakpower_norm(lmi) = TFR_ST(peakF(lmi),peakT(lmi));
  512. end
  514. evBndsF(reg_i,:) = mean(Ffwhm,1);
  515. evBndsT(reg_i,:) = mean(Tfwhm,1);
  516. evPeakF_inds(reg_i) = round(mean(peakF));
  517. evPeakT_inds(reg_i) = round(mean(peakT));
  518. evPeakpower_norm(reg_i) = mean(peakpower_norm);
  519. end
  521. % 12 column matrix with 1. trial index, 2. trial class, 3. maxima frequency, 4. lowerbound frequency, 5. upperbound frequency, 6. frequency span, ...
  522. % 7. maxima timing, 8. event onset timing, 9. event offset timing, 10. event duration, 11. maxima power, 12. maxima/median power, ...
  523. spectralEvents = [spectralEvents; ti*ones(size(evPeakpower)) classLabels(ti)*ones(size(evPeakpower))...
  524. fVec(evPeakF_inds)' evBndsF tVec(evPeakT_inds)' evBndsT evPeakpower evPeakpower_norm];
  525. end
  526. end
  528. end

spectralevents_find.m at commit cd0c83d, under BSD-3-Clause · at the source

Overview

Authors: Christopher J. Black1, Carl Y. Saab2, David A. Borton1,3,4
ORCID iDs: David A. Borton
  1. School of Engineering, Brown University, Providence, RI 02912, USA
  2. Biomedical Engineering, Lerner Research Institute, Cleveland Clinic, Cleveland, OH 44195, USA
  3. Carney Institute for Brain Science, Brown University, Providence, RI 02912, USA
  4. Center for Neurorestoration and Neurotechnology, Rehabilitation R&D Service, Department of Veterans Affairs, Providence, RI 02908, USA
Institutions: Brown University (United States); Cleveland Clinic (United States); United States Department of Veterans Affairs (United States); Providence VA Medical Center (United States); Rehabilitation Research and Development Service (United States)
Journal: iScience, volume 29, issue 9, article 117355
Dates: received 10 October 2025; accepted 11 August 2026; published online 26 August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.isci.2026.117355 · PMID 42699307 · PMCID PMC13544319 · OpenAlex W4361814259
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2023.03.30.534945
Status: code verified
Categories: systems (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Spectral & time-frequency
Keywords: somatosensory cortex, gamma-band activity, sensory modality, sensory discrimination
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIH (MH115895); NINDS; BRAIN Initiative (1R01NS108414-01)
Citations: not cited yet (Europe PMC); 46 references in the paper
Research resources: Adult male mice RRID:IMSR_JAX:017769, Adult male mice RRID:IMSR_JAX:024109

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 3 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, “Spectral event 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;
  • 3 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:

  • it points to a dataset: OSF gu3jr
  • it says that the data are available on request
  • it says that the code is available on request

Read it in the paper: doi.org/10.1016/j.isci.2026.117355.

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

  • Authors: added David A. Borton (0000-0003-0710-3005); removed David A. Borton

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 4 keywords, 3 funders, 46 references, 2 RRIDs.

Cite

This paper

Black, C. J., Saab, C. Y., & Borton, D. A. (2026). Transient gamma events delineate somatosensory modality in S1. iScience, 29(9), 117355. https://doi.org/10.1016/j.isci.2026.117355

BibTeX

@article{black2026transient,
author = {Black, Christopher J. and Saab, Carl Y. and Borton, David A.},
title = {{Transient gamma events delineate somatosensory modality in S1}},
journal = {iScience},
year = {2026},
month = aug,
volume = {29},
number = {9},
pages = {117355},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.117355},
url = {https://doi.org/10.1016/j.isci.2026.117355},
pmid = {42699307},
pmcid = {PMC13544319}
}

RIS

TY - JOUR
AU - Black, Christopher J.
AU - Saab, Carl Y.
AU - Borton, David A.
TI - Transient gamma events delineate somatosensory modality in S1
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/08/26
VL - 29
IS - 9
SP - 117355
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.117355
UR - https://doi.org/10.1016/j.isci.2026.117355
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.117355",
"type": "article-journal",
"title": "Transient gamma events delineate somatosensory modality in S1",
"container-title": "iScience",
"author": [
{
"family": "Black",
"given": "Christopher J."
},
{
"family": "Saab",
"given": "Carl Y."
},
{
"family": "Borton",
"given": "David A."
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "9",
"page": "117355",
"DOI": "10.1016/j.isci.2026.117355",
"PMID": "42699307",
"PMCID": "PMC13544319",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.117355",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
26
]
]
}
}

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.1097/j.pain.0000000000004044 [code]
No effect of rhythmic visual stimulation on experimental pain perception.
Journal: Pain
In common: SciPy, Matplotlib, NumPy, 4 references
[2] doi:10.1016/j.isci.2026.116152 [code]
Integrating multidimensional nociceptive-related cortical features for unsupervised assessment of anesthesia states in rats.
Journal: iScience
In common: 5 references
[3] doi:10.3390/bioengineering13070793
Vibrotactile Stimulation Encoded by Beta and Gamma Bands Varies with Locations on the Upper Limbs.
Journal: Bioengineering (Basel, Switzerland)
In common: 3 references
[4] doi:10.1371/journal.pcbi.1014378 [code]
A mean-field model of neural networks with PV and SOM interneurons reveals connectivity-based mechanisms of gamma oscillations.
Journal: PLoS computational biology
In common: SciPy, Matplotlib, NumPy, 2 references
[5] doi:10.1093/braincomms/fcag145 [code]
GABA&lt;sub&gt;A&lt;/sub&gt; binding correlates with high-frequency EEG: a possible proxy for depolarization in traumatic brain injury.
Journal: Brain communications
In common: SciPy, Matplotlib, NumPy, 2 references
[6] doi:10.1016/j.celrep.2026.117845 [code]
Hippocampal skill memory expansion drives online performance dynamics during skill learning.
Journal: Cell reports
In common: Image Processing Toolbox, SciPy, Matplotlib, 1 other tool, 1 reference
[7] 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: Image Processing Toolbox, SciPy, Matplotlib, 1 other tool, 1 reference
[8] doi:10.1371/journal.pbio.3003824 [code]
Flexible goal learning involves coordinated population activity in dCA1 and medial orbitofrontal cortex.
Journal: PLoS biology
In common: SciPy, Matplotlib, NumPy, systems, 1 reference
[9] doi:10.7554/elife.108408 [code]
Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: Image Processing Toolbox, SciPy, Matplotlib, 1 other tool, systems
[10] doi:10.1016/j.isci.2026.117375 [code]
Motor priming is associated with widespread recruitment into neural ensembles and more rapid ensemble transitions.
Journal: iScience
In common: Image Processing Toolbox, SciPy, Matplotlib, 1 other tool, systems

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.