OSCR

Quantifying Correlogram Shape to Analyze Neuronal Firing Dynamics Recorded in TBI-on-a-Chip.

Code ↔ Paper

12 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 12 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Correlogram Shape Quantification › Quantifying and Classifying Firing Pattern ↔ Analysis Execution/summarizeCorrelograms.m, lines 299–382 · score 0.97 · beta band, gamma band, theta band, delta band, 13–30 Hz, 8–12 Hz
  2. [2] § Methods › Correlogram Shape Quantification › Quantifying and Classifying Firing Pattern ↔ Analysis Execution/calculateCorrelogramMetrics.m, lines 1–147 · score 0.77 · uniform distribution, comparison signal, findpeaks, reference signal, loess, prominence
  3. [3] § Methods › Correlogram Shape Quantification › Quantifying and Classifying Signal Dependence ↔ Analysis Execution/calculateCorrelogramMetrics.m, lines 1–147 · score 0.74 · chi2cdf, chi squared, uniform distribution, cutoff, threshold, nonuniform
  4. [4] § Methods › Correlogram Shape Quantification › Quantifying and Classifying Firing Pattern ↔ Analysis Execution/plotCorrelograms.m, the whole file · a weak match · score 0.69 · smoothed correlogram, comparison signal, reference signal, loess, smoothdata, window
  5. [5] § Methods › Correlogram Shape Quantification ↔ Analysis Execution/summarizeCorrelograms.m, lines 893–987 · score 0.68 · follower strength classifications, nonuniform correlograms, uniformity classifications, Empty, leader, peak
  6. [6] § Methods › Algorithm Pipeline ↔ Analysis Execution/removeBadSignals.m, the whole file · a weak match · score 0.67 · noisy signals, eliminate signals, prompts, empty, rasters, plx
  7. [7] § Results › Correlogram Classification Heat Maps › Bicuculline ↔ Analysis Execution/summarizeCorrelograms.m, lines 204–238 · score 0.64 · fairly weak, fairly strong, strong leader, week, intermediate, follower
  8. [8] § Results › Summarizing Distributions by Classifying Correlogram Metrics ↔ Analysis Execution/summarizeCorrelograms.m, lines 204–238 · score 0.63 · follower nature, correlogram lead, strong leader, weak, fraction, classifies
  9. [9] § Methods › MEA Recording to Validate Algorithm Outputs ↔ Main_AnalyzeRecordingData.m, lines 16–24 · score 0.63 · bicuculline methiodide, impact injuries, permission, reused, publication, network
  10. [10] § Results › Violin Plots to Visualize Correlogram Metric Distributions ↔ Analysis Execution/Violin.m, lines 1–112 · score 0.62 · kernel density, Bastian, Bechtold, Hintze, Nelson, violinplot
  11. [11] § Results › Tracking Changes in Correlogram Classification › Alkalosis ↔ Analysis Execution/summarizeCorrelograms.m, lines 893–987 · score 0.61 · follower strength classification, nonuniform correlograms, correlogram classifications, leader, peaks
  12. [12] § Methods › Correlogram Shape Quantification › Quantifying and Classifying Firing Pattern ↔ Main_AnalyzeRecordingData.m, lines 59–75 · score 0.54 · bin widths, peak locations, threshold, Correlogram peak, classified, signals

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 · 987 lines · 61 KB · MIT · 5 matches

  1. function summarizeCorrelograms(sw, plotProps, AnalysisRegions, InjuryIndicies, ProcessedData, RecordingMetrics, Correlograms)
  2. %This function classifies correlograms based on metrics (uniformity, peak
  3. %count, area left of zero) and plots classification-related outputs
  5. %% Find where recording regions occurr relative to injuries or treatments
  6. if InjuryIndicies.NumberOfInjuriesOrTreatments > 0 %if there was an injury/treatment in the history of the recording
  8. InjuryInds = zeros(3,InjuryIndicies.NumberOfInjuriesOrTreatments); %initialize array to save x coordinate of injuries/treatments for the violin plots (row 1) as well as the shading width (row 2) to indicate where the region occurred, and flag whether the injury occurred before the start of the recording (row 3)
  10. for j = 1: InjuryIndicies.NumberOfInjuriesOrTreatments %cycle through treatments/injuries
  12. relativeStartTimez = AnalysisRegions.Minutes-InjuryIndicies.InjuryStartMin(j); %get the time difference between each region and the injury start
  13. relativeEndTimez = AnalysisRegions.Minutes-InjuryIndicies.InjuryEndMin(j); %get the time difference between each region and the injury end
  15. endNegs = relativeEndTimez<0; %find where a region started (col 1) or ended (col 2) before an injury/treatment end
  16. endPos = relativeEndTimez>0; %find where a region started (col 1) or ended (col 2) after an injury/treatment end
  17. startNegs = relativeStartTimez<0; %find where a region started (col 1) or ended (col 2) before an injury/treatment end
  19. startNegSum = sum(startNegs,2); %if this sum is 2, then the region started and ended before the injury/treatment started.
  20. endNegSum = sum(endNegs,2); %if this sum is 2, then the region started and ended before the injury/treatment ended.
  22. endPosSum = sum(endPos,2); %if this sum is 2, then the region started and ended after the injury/treatment ended.
  24. simultaneousRegAndInj = endNegSum ~= startNegSum; %find whether any region was both in and out of the treatment range/simultaneous with the treatment
  26. if any(simultaneousRegAndInj) %if any row had both positive and negative times/occurred during an injury
  27. injuryIndj = find(simultaneousRegAndInj);
  28. shadeWidthj = 0.4;
  29. startFlagj = 0; %flag=1 if the injury occurred before the recording start
  31. else %otherwise the regions are either completely after or completely after the given injury/treatment
  32. locStartBetween = find(endNegSum==0,1,'first'); %find the first region after the injury/treatment
  33. locEndBetween = find(endPosSum==0,1,'last'); %find the last region before the injury/treatment
  35. if ~isempty(locEndBetween) && ~isempty(locStartBetween) %if the injury/treatment was during the recording
  36. injuryIndj = locStartBetween-(locStartBetween - locEndBetween)/2; %the index of the injury/treatment is between these two regions
  37. shadeWidthj = 0.03;
  38. startFlagj = 0; %flag=1 if the injury occurred before the recording start
  40. elseif isempty(locEndBetween) %if the treatment/injury was before the recording started or after all selected regions
  41. if InjuryIndicies.InjuryStartMin(j) > 0 %injury was after all analyzed regions
  42. injuryIndj = locStartBetween-.5; %the index of the injury/treatment is between these two regions
  43. shadeWidthj = 0.03;
  44. startFlagj = 0; %flag=1 if the injury occurred before the recording startstartFlagj = 0; %flag=1 if the injury occurred before the recording start
  45. else %otherwise injury or treatment before all regions
  46. injuryIndj = locStartBetween/2; %the index of the injury/treatment is between these two regions
  47. if injuryIndj == 0.5; injuryIndj = 0.3; end %shift further from violins if injury/treatment before recording start
  48. shadeWidthj = 0.03;
  49. startFlagj = 1; %flag that the injury occurred before the recording start
  50. end
  51. elseif isempty(locStartBetween) %if the treatment/injury was after the region
  52. injuryIndj = locEndBetween+ 0.5; %the index of the injury/treatment is between these two regions
  53. shadeWidthj = 0.03;
  54. startFlagj = 0; %flag=1 if the injury occurred before the recording start
  55. end
  57. end %end check for simultaneous regions and injuries/treatments
  59. if length(injuryIndj)==1 %if there is only one region, or one space regions before which injury occurred..
  60. InjuryInds(1,j) = injuryIndj; %save the position of injury/treatment j in the list
  61. InjuryInds(2,j) = shadeWidthj; %save the width of the shading box for injury/treatment j in the list
  62. InjuryInds(3,j) = startFlagj; %save the width of the shading box for injury/treatment j in the list
  64. else %otherwise, the injury was administered over multiple regions...
  65. InjuryInds(1,j) = (injuryIndj(end)-injuryIndj(1))/2+injuryIndj(1); %save the position of injury/treatment j in the list
  66. numRegionsInvolved = length(injuryIndj); %get the number of regions in the treatment (to determine shading width)
  67. InjuryInds(2,j) = numRegionsInvolved/2; %save the width of the shading box for injury/treatment j in the list
  68. InjuryInds(3,j) = startFlagj; %save the width of the shading box for injury/treatment j in the list
  69. end
  71. end %end cycle through treatments/injuries
  73. end %end check if there was an injury/treatment
  75. %% Get recording zones
  76. %Zones = times relative to injury or treatment (i.e. before treatment, between treatments, after treatment)
  77. groupOrder = ProcessedData.RegionStats.groupOrder;
  78. timeRange = ProcessedData.RegionStats.timeRange;
  80. % Sort the analysis regions into each recording zone
  81. Zonez = zeros(1,AnalysisRegions.numRegions); %row = zone, column = region
  82. for zz = 1:size(timeRange,2) %cycle through zones...
  84. startTimeLimit = timeRange(1,zz); %the zone starts mid recording, so just look for anything from zone start to zone end
  85. endTimeLimit = timeRange(2,zz); %the zone ends mid recording as well
  87. %Find analysis regions within the zone
  88. %Note that any regions that are part of multiple zones are not marked as being part of either zone
  89. %however, such regions will be used as comparisons against those in the zone later
  91. afterZoneStart = AnalysisRegions.Seconds>=startTimeLimit; %find analysis region times after zone start
  92. afterZoneStart = sum(afterZoneStart,2)==2; %get analysis regions where the region start and end are both after the start of the zone
  94. beforeZoneEnd = AnalysisRegions.Seconds<=endTimeLimit; %find analysis region times before zone end
  95. beforeZoneEnd = sum(beforeZoneEnd,2)==2; %get analysis regions where the region start and end are both before the end of the zone
  97. inZoneInds = afterZoneStart&beforeZoneEnd; %find regions contained within the zone
  99. if ~any(inZoneInds==1) %if the zone is smaller than a single analysis region, ensure that those regions (where the treatmment is administered) are added to the zone
  100. disp([' Zone ' num2str(zz) ' (' char(groupOrder{zz}) ') began and ended in a single region or in the middle of two adjacent regions.'])
  102. afterZoneStart = AnalysisRegions.Seconds>=startTimeLimit; %find analysis region times after zone start
  103. afterZoneStart = sum(afterZoneStart,2)==1; %get analysis regions where the region start and end are both after the start of the zone
  105. beforeZoneEnd = AnalysisRegions.Seconds<=endTimeLimit; %find analysis region times before zone end
  106. beforeZoneEnd = sum(beforeZoneEnd,2)==1; %get analysis regions where the region start and end are both before the end of the zone
  108. zoneLimitz = find(afterZoneStart):find(beforeZoneEnd);
  109. inZoneInds = false(size(afterZoneStart)); inZoneInds(zoneLimitz) = true; %find regions contained within the zone
  110. end
  111. Zonez(inZoneInds) = zz; %mark which zone the region belongs to
  113. end %end cycle through zones
  115. numZones = length(unique(Zonez(Zonez>0))); %count the number of recording zones
  119. %% Get correlogram data arrays
  121. %Determine whether to exclude autocorrelograms from the analysis based on user preferences
  122. if plotProps.includeAutocorrelogramsInStatistics == 1 %if user wants to include autocorrelograms in the violins...
  123. b1 = zeros(size(RecordingMetrics.CorrelogramMetrics.Region1.leaderProbMat)); %include all indicies
  124. else %otherwise, only include crosscorrelograms in violins
  125. %get the autocorrelograms (same number as the signal count, unless only select references are specified)
  126. b1=eye(size(RecordingMetrics.CorrelogramMetrics.Region1.leaderProbMat)); %find the matrix diagonal/autocorrelograms
  127. end
  129. %If only selecting a subset of signals as the reference signal, take only the relevent indicies (of those comparisons only)
  130. b2=ones(ProcessedData.CellCount,ProcessedData.CellCount); b2(:,RecordingMetrics.CorrelogramMetrics.refIndicies)=0;
  132. b=(b1 & ~b2); autoCorNum = sum(b(:)); %count how many correlograms (if any) are being excluded
  134. %Initialize arrays
  135. leaderProbArray = zeros(ProcessedData.CellCount.*length(RecordingMetrics.CorrelogramMetrics.refIndicies)-autoCorNum,AnalysisRegions.numRegions);
  136. numPeaksArray = zeros(ProcessedData.CellCount.*length(RecordingMetrics.CorrelogramMetrics.refIndicies)-autoCorNum,AnalysisRegions.numRegions);
  137. uniformityArray = zeros(ProcessedData.CellCount.*length(RecordingMetrics.CorrelogramMetrics.refIndicies)-autoCorNum,AnalysisRegions.numRegions);
  138. uniformitypArray = zeros(ProcessedData.CellCount.*length(RecordingMetrics.CorrelogramMetrics.refIndicies)-autoCorNum,AnalysisRegions.numRegions);
  140. %For peak positions (which provide frequencies for any rhythmic activity), each region can have a different number of peaks. Therefore, regions with less than the maximum number of peaks need to be padded with dummy numbers so that they are all the same size for violin plot
  141. MaxNumPk = 0;
  142. for rr = 1:AnalysisRegions.numRegions %find the maximum possible number of peaks for all regions
  143. reg_field = strcat('Region',num2str(rr)); %convert the name of the region to a variable name to call the correct correlogram substructure
  144. MaxNumPk = max(MaxNumPk,max(RecordingMetrics.CorrelogramMetrics.(reg_field).numberOfCorrelogramPeaks));
  145. end
  146. PkArray = NaN(MaxNumPk,AnalysisRegions.numRegions); %initialize peak frequency array
  148. SourceCorrelogram = zeros(size(PkArray)); %initialize array to store an ID for the correlogram in which each peak is located
  150. %% Populate the correlogram arrays for each analysis region
  151. for rr = 1:AnalysisRegions.numRegions %for each region...
  153. reg_field = strcat('Region',num2str(rr)); %get the region name to access the data structure
  154. if plotProps.includeAutocorrelogramsInStatistics == 1 %if user wants to include autocorrelograms in the violins
  155. leaderProbArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).leaderProb'; %store the metric in the data structure
  156. numPeaksArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).numberOfCorrelogramPeaks'; %store the metric in the data structure
  157. uniformityArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).uniformityTest'; %store the metric in the data structure
  158. b1 = zeros(size(RecordingMetrics.CorrelogramMetrics.(reg_field).leaderProbMat));
  159. else
  160. b1=eye(size(RecordingMetrics.CorrelogramMetrics.(reg_field).leaderProbMat)); %find the matrix diagonal/autocorrelograms
  161. end
  163. b2=ones(ProcessedData.CellCount,ProcessedData.CellCount); b2(:,RecordingMetrics.CorrelogramMetrics.refIndicies)=0;
  164. b = b1 | b2; %get only the reference signals you want (if specific references were selected);
  166. %get all cross-correlogram properties (not including autocorrelograms or excluded references if specified)
  167. leaderProbArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).leaderProbMat(~b); %store the metric in the data structure
  168. numPeaksArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).numberOfCorrelogramPeaks(~b); %store the metric in the data structure
  169. uniformityArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).uniformityTest(~b); %store the metric in the data structure
  170. uniformitypArray(:,rr) = RecordingMetrics.CorrelogramMetrics.(reg_field).uniformityTestPValue(~b); %store the metric in the data structure
  172. if plotProps.classifyPeakTimesAsFrequencies == 1 %if classifying peak times according to frequency (Hz, f = 1/t)...
  173. PeakLocationArray = cell2mat(cellfun(@(x) 1./abs(x), RecordingMetrics.CorrelogramMetrics.(reg_field).correlogramPeakLocations(~b),'UniformOutput',false))'; %list peak frequencies
  174. PeakAtZero = isinf(PeakLocationArray); PeakLocationArray = PeakLocationArray(~PeakAtZero); %remove peaks that occurred at t=0, which have f = 1/0 = infinity (which is not physiological)
  175. PkArray(1:length(PeakLocationArray),rr) = PeakLocationArray; %store the metric in the data structure
  176. else %otherwise classify peak times accorsing to time (s)
  177. PeakLocationArray = cell2mat(cellfun(@(x) abs(x), RecordingMetrics.CorrelogramMetrics.(reg_field).correlogramPeakLocations(~b),'UniformOutput',false))'; %list peak frequencies
  178. PkArray(1:length(PeakLocationArray),rr) = PeakLocationArray; %store the metric in the data structure
  179. end
  181. %Store the ID of the correlogram from which each peak originated
  182. startInd = 1; %start index for source correlogram IDs (since the number of peaks can be > 1 for any correlogram, the peak frequency array is larger than the correlogram arrays. Therefore, it is necessary to track which correlogram contained each peak)
  183. for i = 1:ProcessedData.CellCount^2 %cycle through correlograms...
  184. N = RecordingMetrics.CorrelogramMetrics.(reg_field).numberOfCorrelogramPeaks(i); %Get the number of correlogram peaks
  185. N2 = sum(isinf(1./(abs(RecordingMetrics.CorrelogramMetrics.(reg_field).correlogramPeakLocations{i})))); %determine if f = Inf (unphysiological) for any of the peaks
  187. if isnan(N) || N > N2 %if there is at least one peak not equal to infinity...
  188. if isnan(N) %set the value of N if the correlogram is empty
  189. N = 1; N2 = 0;
  190. end
  192. SourceCorrelogram(startInd:startInd+(N-N2)-1,rr) = i; %save correlogram source ID in the matrix
  193. startInd = startInd+(N-N2); %advance to the start for the next correlogram ID
  195. end %end check that there is at least one peak not equal to infinity
  197. end %end cycle through correlograms
  199. end %end cycle through regions
  201. PkArray(PkArray==0)=NaN;
  204. %% Classify leader/follower nature of the correlograms
  205. leaderProbGroupBoundaries = 0:.2:1;
  206. leaderClassificationNames = {'Weak', 'Fairly Weak', 'Intermediate', 'Fairly Strong', 'Strong'};
  207. NumberOfleaderClassificationGroups = length(leaderClassificationNames); %get all possible groups for nonempty correlograms
  209. fracInProbGroup = NaN(size(leaderProbGroupBoundaries,2)-1,AnalysisRegions.numRegions); %Row = grouping, Col = Recording Region
  210. LFGroupAssignments = NaN(size(leaderProbArray)); %initialize array to store leader/follower classifications
  212. StrongWeakClass = abs(leaderProbArray - 0.5)./0.5; %Bin into five categories of weak->strong leader/follower characteristic. A value of 1 is strong, 0 is weak
  214. for g = 1:length(leaderProbGroupBoundaries)-1
  215. %Get bounds on the probabilities in range
  216. a = leaderProbGroupBoundaries(g); b = leaderProbGroupBoundaries(g+1);
  218. %Find the values of how strong/week a given correlogram leads/follows within in the given range (a - b)
  219. if g == 1 %at the start of category listing, so pad the zero (lower group bound) to avoid machine/rounding error
  220. ProbInBin = (StrongWeakClass>=a-1 & StrongWeakClass<b); %find leader/follower strength classifications within the given range
  221. elseif b~=1 %as long as b is not the end of the list
  222. ProbInBin = (StrongWeakClass>=a & StrongWeakClass<b); %find leader/follower strength classifications within the given range
  223. else %otherwise, end of the category listing, so be sure to include the highest limit (don't reserve for the next group as with non-end elements of the categories)
  224. %at the start of category listing, so pad the 1 (upper group bound) to avoid machine/rounding error
  225. ProbInBin = (StrongWeakClass>=a & StrongWeakClass<=b+1); %find leader/follower strength classifications within the given range
  226. end
  228. LFGroupAssignments(ProbInBin) = g; %save the logical array marking which correlograms fall into which bin
  230. fracInProbGroup(g,1:end) = sum(ProbInBin,1); %sum up how many correlograms fit in the given bin
  232. end
  234. %Do not count correlograms with 0 events in the total
  235. EmptyCorrelograms = isnan(leaderProbArray); %Find correlograms with 0 events
  236. totalNumberOfCorrelograms = sum(~EmptyCorrelograms,1); %Get the actual/non-empty count
  238. fracInProbGroup = fracInProbGroup./totalNumberOfCorrelograms; %Get the fraction of correlograms in each bin
  241. %% Classify peak counts of the correlograms
  242. pkCountGroupBoundaries = min(numPeaksArray(:)):1:min(plotProps.maxPeakCountBeforeNoise,max(numPeaksArray(:))); %bin peak counts from the minimum to a user-specified maximum or the true maximum (whichever max value is smaller)
  243. pkCountNames = string(pkCountGroupBoundaries); %group names are just the peak count
  245. if plotProps.maxPeakCountBeforeNoise < max(numPeaksArray(:)) %the last group/peak count is a user-specified maximum, which may be smaller than the true maximum, so add "or more" to the name if this is the case
  246. pkCountNames{end} = strcat(pkCountNames{end},' or More');
  247. end
  248. NumberOfpkCountClassificationGroups = length(pkCountNames); %get all possible groups for nonempty correlograms
  250. fracInPeakGroup = NaN(size(pkCountGroupBoundaries,2)-1,AnalysisRegions.numRegions); %Row = grouping, Col = Recording Region
  251. PCGroupAssignments = NaN(size(numPeaksArray)); %initialize array to store peak count classifications
  253. for g = 1:length(pkCountGroupBoundaries)
  254. %Get bounds on the probabilities in range
  255. if g == length(pkCountGroupBoundaries)
  256. ProbInBin = numPeaksArray >= pkCountGroupBoundaries(g); %find correlograms with the given number of peaks or more if at the end of the group list
  257. else
  258. ProbInBin = numPeaksArray == pkCountGroupBoundaries(g); %find correlograms with the given number of peaks
  259. end
  261. PCGroupAssignments(ProbInBin) = g; %save the logical array marking which correlograms fall into which bin
  263. fracInPeakGroup(g,1:end) = sum(ProbInBin,1); %sum up how many correlograms fit in the given bin
  265. end
  267. %Do not count correlograms with 0 events in the total
  268. EmptyCorrelograms = isnan(numPeaksArray); %Find correlograms with 0 events
  269. totalNumberOfCorrelograms = sum(~EmptyCorrelograms,1); %Get the actual/non-empty count
  271. fracInPeakGroup = fracInPeakGroup./totalNumberOfCorrelograms; %Get the fraction of correlograms in each bin
  273. %% Uniformity is already classified as uniform and nonuniform, so just get the fraction of uniform and nonuniform in each data analysis region
  274. unifGroupBoundaries = [0 1];
  275. unifNames = {'Nonuniform','Uniform'};
  276. NumberOfunifClassificationGroups = length(unifGroupBoundaries); %get all possible groups for nonempty correlograms
  278. fracInUnifGroup = NaN(size(unifGroupBoundaries,2)-1,AnalysisRegions.numRegions); %Row = grouping, Col = Recording Region
  279. unifGroupAssignments = NaN(size(uniformityArray)); %initialize array to store peak count classifications
  281. for g = 1:length(unifGroupBoundaries)
  282. %Get bounds on the probabilities in range
  283. ProbInBin = uniformityArray == unifGroupBoundaries(g); %find correlograms with the given uniformity classification
  285. unifGroupAssignments(ProbInBin) = g; %save the logical array marking which correlograms fall into which bin
  287. fracInUnifGroup(g,1:end) = sum(ProbInBin,1); %sum up how many correlograms fit in the given bin
  289. end
  291. %Do not count correlograms with 0 events in the total
  292. EmptyCorrelograms = isnan(uniformityArray); %Find correlograms with 0 events
  293. totalNumberOfCorrelograms = sum(~EmptyCorrelograms,1); %Get the actual/non-empty count
  295. fracInUnifGroup = fracInUnifGroup./totalNumberOfCorrelograms; %Get the fraction of correlograms in each bin
  299. %% Classify correlogram peak times or peak frequencies
  300. % EEG Categories (if classifying by frequency)
  301. % Delta: f < 4
  302. % Theta: 4 <= f <= 7
  303. % Alpha: 8 <= f <= 12
  304. % Beta: 13 <= f <= 30
  305. % Gamma: f > 32
  307. if plotProps.classifyPeakTimesAsFrequencies == 1 %if classifying peak times according to frequency (Hz, f = 1/t)...
  308. startF = [0 4 7 8 12 13 30 32]; %list the minimum frequency of each category
  309. endF = [4 7 8 12 13 30 32 max(PkArray(:))+10]; %list the maximum frequency of each category
  310. freqClassificationNames = {'Delta Band (<4 Hz)', 'Theta Band (4-7 Hz)', 'Between Theta and Alpha (7-8 Hz)', 'Alpha Band (8-12 Hz)', 'Between Alpha and Beta (12-13 Hz)', 'Beta Band (13-30 Hz)', 'Between Beta and Gamma (30-32 Hz)', 'Gamma Band (>32 Hz)'};
  311. NumberOfFreqClassificationGroups = length(freqClassificationNames); %count the number of frequency categories
  312. else
  313. stpowr = [-Inf floor(log10(Correlograms.CorrelogramBins(end)-Correlograms.CorrelogramBins(end-1))):1:ceil(log10(max(PkArray(:))))-1]; %get the log10 power of time bin starts
  314. startF = 10.^stpowr; %list the minimum time of each category
  315. endpowr = [floor(log10(Correlograms.CorrelogramBins(end)-Correlograms.CorrelogramBins(end-1))):1:ceil(log10(max(PkArray(:))))]; %get the log10 power of time bin ends
  316. endF = 10.^endpowr; %list the maximum time of each category
  318. freqClassificationNames = repmat({' '},1,length(stpowr)); %initialize array to save time grouping names
  319. for i = 1:length(stpowr)
  320. if i == 1; ind = i+1; else; ind = i; end %get index to access the start time name
  322. %get the start and end time strings
  323. if stpowr(ind) == -12; sttstr = '1 picosecond'; elseif stpowr(ind) == -11; sttstr = '10 picoseconds'; elseif stpowr(ind) == -10; sttstr = '100 picoseconds'; elseif stpowr(ind) == -9; sttstr = '1 nanosecond'; elseif stpowr(ind) == -8; sttstr = '10 nanoseconds'; elseif stpowr(ind) == -7; sttstr = '100 nanoseconds'; elseif stpowr(ind) == -6; sttstr = '1 microsecond'; elseif stpowr(ind) == -5; sttstr = '10 microseconds'; elseif stpowr(ind) == -4; sttstr = '100 microseconds'; elseif stpowr(ind) == -3; sttstr = '1 millisecond'; elseif stpowr(ind) == -2; sttstr = '10 milliseconds'; elseif stpowr(ind) == -1; sttstr = '100 milliseconds'; elseif stpowr(ind) == 0; sttstr = '1 second'; elseif stpowr(ind) == 1; sttstr = '10 seconds'; elseif stpowr(ind) == 2; sttstr = '100 seconds'; elseif stpowr(ind) == 3; sttstr = '1,000 seconds'; elseif stpowr(ind) == 4; sttstr = '10,000 seconds'; elseif stpowr(ind) == 5; sttstr = '100,000 seconds'; elseif stpowr(ind) == 6; sttstr = '1,000,000 seconds'; end
  324. if endpowr(i) == -12; endtstr = '1 picosecond'; elseif endpowr(i) == -11; endtstr = '10 picoseconds'; elseif endpowr(i) == -10; endtstr = '100 picoseconds'; elseif endpowr(i) == -9; endtstr = '1 nanosecond'; elseif endpowr(i) == -8; endtstr = '10 nanoseconds'; elseif endpowr(i) == -7; endtstr = '100 nanoseconds'; elseif endpowr(i) == -6; endtstr = '1 microsecond'; elseif endpowr(i) == -5; endtstr = '10 microseconds'; elseif endpowr(i) == -4; endtstr = '100 microseconds'; elseif endpowr(i) == -3; endtstr = '1 millisecond'; elseif endpowr(i) == -2; endtstr = '10 milliseconds'; elseif endpowr(i) == -1; endtstr = '100 milliseconds'; elseif endpowr(i) == 0; endtstr = '1 second'; elseif endpowr(i) == 1; endtstr = '10 seconds'; elseif endpowr(i) == 2; endtstr = '100 seconds'; elseif endpowr(i) == 3; endtstr = '1,000 seconds'; elseif endpowr(i) == 4; endtstr = '10,000 seconds'; elseif endpowr(i) == 5; endtstr = '100,000 seconds'; elseif endpowr(i) == 6; endtstr = '1,000,000 seconds'; end
  326. if i == 1 %get the first entry in the categories
  327. tstr = [sttstr ' or less'];
  328. else %now get the non-start time ranges
  329. startchar = strfind(sttstr,' '); startchar = startchar(1); %get the index when the first number stops
  330. if startchar == 4 %if the first number is 100, use the full start name
  331. tstr = [sttstr ' - ' endtstr];
  332. else %otherwise, you only need to write the units once
  333. tstr = [sttstr(1:startchar-1) ' - ' endtstr];
  334. end
  335. end %end get the category range
  336. freqClassificationNames{i} = tstr; %save the category in the time classification name array
  337. end
  339. NumberOfFreqClassificationGroups = length(freqClassificationNames); %count the number of peak time categories
  341. end
  343. %Initialize peak frequency (or time) classification arrays
  344. fracInFreqGroup = NaN(NumberOfFreqClassificationGroups,AnalysisRegions.numRegions); %Row = fraction of correlogram peaks in grouping, Col = Recording Region
  345. meanFreqInFreqGroup = NaN(NumberOfFreqClassificationGroups,AnalysisRegions.numRegions); %Row = mean frequency within the given group, Col = Recording Region
  346. sdFreqInFreqGroup = NaN(NumberOfFreqClassificationGroups,AnalysisRegions.numRegions); %Row = standard deviation in the frequency within the given group, Col = Recording Region
  347. seFreqInFreqGroup = NaN(NumberOfFreqClassificationGroups,AnalysisRegions.numRegions); %Row = standard deviation in the frequency within the given group, Col = Recording Region
  348. minFreqInFreqGroup = NaN(NumberOfFreqClassificationGroups,AnalysisRegions.numRegions); %Row = min frequency within the given group, Col = Recording Region
  349. maxFreqInFreqGroup = NaN(NumberOfFreqClassificationGroups,AnalysisRegions.numRegions); %Row = max frequency within the given group, Col = Recording Region
  350. FreqGroupAssignments = NaN(size(PkArray)); %initialize array to store peak frequency classifications
  352. for g = 1:NumberOfFreqClassificationGroups %for each classification grouping...
  354. if g == 1
  355. b = PkArray<= endF(g);
  356. else
  357. b = PkArray>startF(g) & PkArray <= endF(g); %find frequencies that fall within the given classification/grouping
  358. end
  359. FreqGroupAssignments(b) = g; %save the logical array marking peaks in the given classification
  361. fracInFreqGroup(g,:) = sum(b,1); %count how many correlogram peaks have a frequency in the given classification
  363. %min and max functions will return an empty array if no peaks are in group g for a given region. This interferes with MATLAB's indexing, so empty regions must be found and ignored for the min and max functions
  364. regionsWithPeaksInGroup = sum(b,1)>0; %make a logical indicating regions with at least one peak in the given classification group
  365. regionsWithGroup = 1:AnalysisRegions.numRegions; regionsWithGroup = regionsWithGroup(regionsWithPeaksInGroup); %list regions with at least one peak in the given classification group
  367. if any(regionsWithPeaksInGroup)
  368. meanFreqInFreqGroup(g,:) = cell2mat(arrayfun(@(col) mean(PkArray(b(:,col),col),'omitnan'),1:AnalysisRegions.numRegions,'UniformOutput',false)); % mean frequency of correlogram peaks in the given classification
  369. medianFreqInFreqGroup(g,:) = cell2mat(arrayfun(@(col) nanmedian(PkArray(b(:,col),col)),1:AnalysisRegions.numRegions,'UniformOutput',false)); % median frequency of correlogram peaks in the given classification
  370. sdFreqInFreqGroup(g,:) = cell2mat(arrayfun(@(col) std(PkArray(b(:,col),col),'omitnan'),1:AnalysisRegions.numRegions,'UniformOutput',false)); % standard deviation in the frequencies of correlogram peaks in the given classification
  371. seFreqInFreqGroup(g,:) = cell2mat(arrayfun(@(col) std(PkArray(b(:,col),col)/sqrt(length(b(:,col))),'omitnan'),1:AnalysisRegions.numRegions,'UniformOutput',false)); % standard error in the frequencies of correlogram peaks in the given classification
  372. minFreqInFreqGroup(g,regionsWithGroup) = cell2mat(arrayfun(@(col) min(PkArray(b(:,col),col),[],'omitnan'),regionsWithGroup,'UniformOutput',false)); % min frequency of correlogram peaks in the given classification
  373. maxFreqInFreqGroup(g,regionsWithGroup) = cell2mat(arrayfun(@(col) max(PkArray(b(:,col),col),[],'omitnan'),regionsWithGroup,'UniformOutput',false)); % max frequency of correlogram peaks in the given classification
  374. end
  376. end
  378. %Do not count correlograms with 0 events in the total for the % of peaks in each classification
  379. EmptyCorrelograms = isnan(PkArray) | PkArray == 0; %Find correlograms with 0 events
  380. totalNumberOfCorrelograms = sum(~EmptyCorrelograms,1); %Get the actual/non-empty count
  382. fracInFreqGroup = fracInFreqGroup./totalNumberOfCorrelograms; %Get the fraction of correlograms in each bin (sum of all rows for each column = 1)
  385. %% Make a structure containing indicies for each classification group
  386. %(for ease of accessing later/to prevent constantly using "find")
  388. for g = 1:NumberOfleaderClassificationGroups
  389. g_name = strcat('G',num2str(g));
  390. FindGroups.LeaderFollower.(g_name) = LFGroupAssignments==g;
  391. end
  392. for g = 1:NumberOfpkCountClassificationGroups
  393. g_name = strcat('G',num2str(g));
  394. FindGroups.PeakCount.(g_name) = PCGroupAssignments==g;
  395. end
  396. for g = 1:NumberOfunifClassificationGroups
  397. g_name = strcat('G',num2str(g));
  398. FindGroups.Unif.(g_name) = unifGroupAssignments==g;
  399. end
  400. for g = 1:NumberOfFreqClassificationGroups
  401. g_name = strcat('G',num2str(g));
  402. FindGroups.FreqClassification.(g_name) = FreqGroupAssignments==g;
  403. end
  405. %% Analyze classification changes between treatments/injuries
  406. %Get the probability that the correlogram now has its current classification, given the classification in the previous recording region
  407. for zz = 1: numZones %cycle through recording zones
  408. zoneName = strcat('Zone',num2str(zz)); %get the zone name to access the data structure
  409. startR = find(Zonez==zz,1,'first'); %get the first region in the zone
  410. endR = find(Zonez==zz,1,'last'); %get the last region in the zone
  412. LFGroupz = zeros(NumberOfleaderClassificationGroups,NumberOfleaderClassificationGroups); %Initialize array to save probability of being in a group given the past group
  413. PCGroupz = zeros(NumberOfpkCountClassificationGroups,NumberOfpkCountClassificationGroups); %Initialize array to save probability of being in a group given the past group
  414. UnifGroupz = zeros(NumberOfunifClassificationGroups,NumberOfunifClassificationGroups); %Initialize array to save probability of being in a group given the past group
  416. %row = previous group, column = probability of current group
  417. for rr = max(2,startR):endR %cycle through regions within the zone...
  419. %Get the probability of a current leader/follower classification given the value of the previous classification
  420. for g = 1:NumberOfleaderClassificationGroups %cycle through past group...
  421. g_name = ['G' num2str(g)]; %get the group name
  422. indArr = FindGroups.LeaderFollower.(g_name); %find correlograms classified in the given group
  424. previouslyInGroup = indArr(:,rr-1); %find correlograms belonging to the given group in the previous analysis region
  425. currentGroup = LFGroupAssignments(previouslyInGroup,rr); %find how those correlograms are classified in the current analysis region
  427. for g2 = 1:NumberOfleaderClassificationGroups %for each current group...
  428. LFGroupz(g,g2) = LFGroupz(g,g2)+sum(currentGroup==g2); %add current groups to histogram counts
  429. end %end cycle through current groups
  431. end %end cycle through past group
  434. %Get the probability of a current peak count given the value of the previous peak count
  435. for g = 1:NumberOfpkCountClassificationGroups %cycle through past group...
  436. g_name = ['G' num2str(g)]; %get the group name
  437. indArr = FindGroups.PeakCount.(g_name); %find correlograms classified in the given group
  439. previouslyInGroup = indArr(:,rr-1); %find correlograms belonging to the given group in the previous analysis region
  440. currentGroup = PCGroupAssignments(previouslyInGroup,rr); %find how those correlograms are classified in the current analysis region
  442. for g2 = 1:NumberOfpkCountClassificationGroups %for each current group...
  443. PCGroupz(g,g2) = PCGroupz(g,g2)+sum(currentGroup==g2); %add current groups to histogram counts
  444. end %end cycle through current groups
  446. end %end cycle through past group
  449. %Get the probability of a current uniformity classification given the value of the previous uniformity classification
  450. for g = 1:NumberOfunifClassificationGroups %cycle through past group...
  451. g_name = ['G' num2str(g)]; %get the group name
  452. indArr = FindGroups.Unif.(g_name); %find correlograms classified in the given group
  454. previouslyInGroup = indArr(:,rr-1); %find correlograms belonging to the given group in the previous analysis region
  455. currentGroup = unifGroupAssignments(previouslyInGroup,rr); %find how those correlograms are classified in the current analysis region
  457. for g2 = 1:NumberOfunifClassificationGroups %for each current group...
  458. UnifGroupz(g,g2) = UnifGroupz(g,g2)+sum(currentGroup==g2); %add current groups to histogram counts
  459. end %end cycle through current groups
  461. end %end cycle through past group
  463. end %end cycle through analysis regions within the zone
  465. ZoneHistograms.LeaderFollower.(zoneName) = LFGroupz./sum(LFGroupz,2); %normalize histogram to probabilities and save in the structure so you can plot histograms for each zone
  466. ZoneHistograms.PeakCount.(zoneName) = PCGroupz./sum(PCGroupz,2); %normalize histogram to probabilities and save in the structure so you can plot histograms for each zone
  467. ZoneHistograms.Unif.(zoneName) = UnifGroupz./sum(UnifGroupz,2); %normalize histogram to probabilities and save in the structure so you can plot histograms for each zone
  469. end %end cycle through recording zones
  472. %% Make Plots of Classification Outputs
  474. %Initialize plot markers for ease of plotting
  475. markerNmArr = {'o','*','x','square','diamond','^','v','>','<','pentagram','hexagram','+','.'}; %Make an array to change the plot marker shapes for each category of classified data
  476. markerSzArr = [plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize*0.4 plotProps.markerSize];
  478. %Make sure there are enough marker indicators for all plots to be made
  479. maxClassificationCount = (1+min(plotProps.maxPeakCountBeforeNoise,max(numPeaksArray(:)))); %The highest number of classifications will be peak counts, so make sure length(markerArray) is larger than the max possible number of peak counts
  480. maxClassificationCount = max(length(ProcessedData.RegionStats.groupOrder),maxClassificationCount); %the highest possible number of lines on that plots will be the number of analysis regions
  482. F = maxClassificationCount/length(markerNmArr); %The highest number of classifications will be peak counts, so make sure length(markerArray) is larger than the max possible number of peak counts
  483. if F > 1 %if the marker array is too short to plot all peak classifications
  484. F = ceil(F); %round up so the array will eventually be longer than needed
  485. for i = 1:F %as many times as needed...
  486. markerNmArr = [markerNmArr markerNmArr]; %lengthen the marker name array as many times as needed
  487. markerSzArr = [markerSzArr markerSzArr];
  488. end
  489. end %end check marker specification is long enough
  493. %Plot the fraction of correlogram peaks in each leader/follower classification
  494. cc = pink(NumberOfleaderClassificationGroups+4); %set plot colors
  496. figure()
  497. %Indicate where injuries/treatments occurred
  498. hold on
  499. if InjuryIndicies.NumberOfInjuriesOrTreatments > 0 %if there was an injury, shade where the injury occurred
  500. for j = 1:InjuryIndicies.NumberOfInjuriesOrTreatments %for all injuries/treatments...
  501. plot([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  502. plot([InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  503. fill([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2 2 -1],[0.6 0.6 0.6],'facealpha',0.2,'LineStyle','none')
  504. if InjuryInds(3,j) ~= 1 %if injury/treatment is mid-recoding
  505. text(InjuryInds(1,j) ,max(fracInProbGroup(:))+.05,InjuryIndicies.InjuryLabels{j},'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  506. else %otherwise injury/treatment is before recording start
  507. text(InjuryInds(1,j) ,max(fracInProbGroup(:))+.05,[InjuryIndicies.InjuryLabels{j} ' Before Recording Start'],'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  508. end
  509. end
  510. end
  511. %Now plot the fraction of correlogram peaks in each classification
  512. h=plot(1:AnalysisRegions.numRegions,fracInProbGroup','.-','LineWidth',plotProps.lineWidth);
  513. hold off
  514. for c = 1:size(h,1)
  515. h(c).Color = cc(c,:);
  516. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  517. end
  518. set(gca,'XTick',1:AnalysisRegions.numRegions)
  519. if isfield(AnalysisRegions,'RegionLabels')
  520. set(gca,'XTickLabel',AnalysisRegions.RegionLabels)
  521. end
  523. xlim([0.5 AnalysisRegions.numRegions+0.5])
  524. xlabel('Recording Region')
  525. ylim([0 1.2*(max(fracInProbGroup(:)))])
  526. ylabel('Fraction of Correlograms with a given Classification')
  527. legend(h,leaderClassificationNames)
  528. set(gca,'FontSize',plotProps.FigureFontSize)
  529. title('Classification of Leader/Follower Nature Through Time')
  530. box on
  531. set(gca,'YTick',0:0.05:1)
  532. ax = gca;
  533. ax.YGrid = 'on';
  537. %Plot the fraction of correlogram peaks in each peak count
  538. cc = pink(NumberOfpkCountClassificationGroups+4); %set plot colors
  540. figure()
  541. %Indicate where injuries/treatments occurred
  542. hold on
  543. if InjuryIndicies.NumberOfInjuriesOrTreatments > 0 %if there was an injury, shade where the injury occurred
  544. for j = 1:InjuryIndicies.NumberOfInjuriesOrTreatments %for all injuries/treatments...
  545. plot([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  546. plot([InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  547. fill([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2 2 -1],[0.6 0.6 0.6],'facealpha',0.2,'LineStyle','none')
  548. if InjuryInds(3,j) ~= 1 %if injury/treatment is mid-recoding
  549. text(InjuryInds(1,j) ,max(fracInPeakGroup(:))+.05,InjuryIndicies.InjuryLabels{j},'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  550. else %otherwise injury/treatment is before recording start
  551. text(InjuryInds(1,j) ,max(fracInPeakGroup(:))+.05,[InjuryIndicies.InjuryLabels{j} ' Before Recording Start'],'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  552. end
  553. end
  554. end
  555. %Now plot the fraction of correlogram peaks in each classification
  556. h=plot(1:AnalysisRegions.numRegions,fracInPeakGroup','.-','LineWidth',plotProps.lineWidth);
  557. hold off
  558. for c = 1:size(h,1)
  559. h(c).Color = cc(c,:);
  560. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  561. end
  562. set(gca,'XTick',1:AnalysisRegions.numRegions)
  563. if isfield(AnalysisRegions,'RegionLabels')
  564. set(gca,'XTickLabel',AnalysisRegions.RegionLabels)
  565. end
  567. xlim([0.5 AnalysisRegions.numRegions+0.5])
  568. xlabel('Recording Region')
  569. ylim([0 1.2*(max(fracInPeakGroup(:)))])
  570. ylabel('Fraction of Correlograms with a given Peak Count')
  571. legend(h,pkCountNames)
  572. set(gca,'FontSize',plotProps.FigureFontSize)
  573. title('Peak Counts Through Time')
  574. box on
  575. set(gca,'YTick',0:0.1:1)
  576. ax = gca;
  577. ax.YGrid = 'on';
  581. %Plot the fraction of correlogram that are uniform and nonuniform
  582. cc = pink(4); %set plot colors
  584. figure()
  585. %Indicate where injuries/treatments occurred
  586. hold on
  587. if InjuryIndicies.NumberOfInjuriesOrTreatments > 0 %if there was an injury, shade where the injury occurred
  588. for j = 1:InjuryIndicies.NumberOfInjuriesOrTreatments %for all injuries/treatments...
  589. plot([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  590. plot([InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  591. fill([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2 2 -1],[0.6 0.6 0.6],'facealpha',0.2,'LineStyle','none')
  592. if InjuryInds(3,j) ~= 1 %if injury/treatment is mid-recoding
  593. text(InjuryInds(1,j) ,max(fracInUnifGroup(:))+.05,InjuryIndicies.InjuryLabels{j},'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  594. else %otherwise injury/treatment is before recording start
  595. text(InjuryInds(1,j) ,max(fracInUnifGroup(:))+.05,[InjuryIndicies.InjuryLabels{j} ' Before Recording Start'],'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  596. end
  597. end
  598. end
  599. %Now plot the fraction of correlogram peaks in each classification
  600. h=plot(1:AnalysisRegions.numRegions,fracInUnifGroup','.-','LineWidth',plotProps.lineWidth);
  601. hold off
  602. for c = 1:size(h,1)
  603. h(c).Color = cc(c,:);
  604. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  605. end
  606. set(gca,'XTick',1:AnalysisRegions.numRegions)
  607. if isfield(AnalysisRegions,'RegionLabels')
  608. set(gca,'XTickLabel',AnalysisRegions.RegionLabels)
  609. end
  611. xlim([0.5 AnalysisRegions.numRegions+0.5])
  612. xlabel('Recording Region')
  613. ylim([0 1.1*(max(fracInUnifGroup(:)))])
  614. ylabel('Fraction of Correlograms with a given Uniformity')
  615. legend(h,unifNames)
  616. set(gca,'FontSize',plotProps.FigureFontSize)
  617. title('Correlogram Uniformity Through Time')
  618. box on
  619. set(gca,'YTick',0:0.1:1)
  620. ax = gca;
  621. ax.YGrid = 'on';
  623. %Plot the fraction of correlogram peaks in each EEG classification
  624. cc = pink(NumberOfFreqClassificationGroups+4); %set plot colors
  626. figure()
  627. %Indicate where injuries/treatments occurred
  628. hold on
  629. if InjuryIndicies.NumberOfInjuriesOrTreatments > 0 %if there was an injury, shade where the injury occurred
  630. for j = 1:InjuryIndicies.NumberOfInjuriesOrTreatments %for all injuries/treatments...
  631. plot([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  632. plot([InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2],'--k','LineWidth',plotProps.lineWidth/2)
  633. fill([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[-1 2 2 -1],[0.6 0.6 0.6],'facealpha',0.2,'LineStyle','none')
  634. if InjuryInds(3,j) ~= 1 %if injury/treatment is mid-recoding
  635. text(InjuryInds(1,j) ,max(fracInFreqGroup(:))+.05,InjuryIndicies.InjuryLabels{j},'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  636. else %otherwise injury/treatment is before recording start
  637. text(InjuryInds(1,j) ,max(fracInFreqGroup(:))+.05,[InjuryIndicies.InjuryLabels{j} ' Before Recording Start'],'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  638. end
  639. end
  640. end
  641. %Now plot the fraction of correlogram peaks in each classification
  642. h=plot(1:AnalysisRegions.numRegions,fracInFreqGroup','.-','LineWidth',plotProps.lineWidth);
  643. hold off
  644. for c = 1:size(h,1)
  645. h(c).Color = cc(c,:);
  646. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  647. end
  648. set(gca,'XTick',1:AnalysisRegions.numRegions)
  649. if isfield(AnalysisRegions,'RegionLabels')
  650. set(gca,'XTickLabel',AnalysisRegions.RegionLabels)
  651. end
  652. xlim([0.5 AnalysisRegions.numRegions+0.5])
  653. xlabel('Recording Region')
  654. ylim([0 1.2*(max(fracInFreqGroup(:)))])
  655. if plotProps.classifyPeakTimesAsFrequencies == 1
  656. ylabel('Fraction of Correlogram Peaks in Each Frequency Band')
  657. title('Peak Frequency Classification Through Time')
  658. else
  659. ylabel('Fraction of Correlogram Peaks in Each Time Range')
  660. title('Peak Location Classification Through Time')
  661. end
  662. legend(h,freqClassificationNames)
  663. set(gca,'FontSize',plotProps.FigureFontSize)
  664. box on
  665. set(gca,'YTick',0:0.1:1)
  666. ax = gca;
  667. ax.YGrid = 'on';
  670. %Plot the average plus/minus standard deviation of the frequency of each band
  671. errorbarx = zeros(AnalysisRegions.numRegions,NumberOfFreqClassificationGroups); %set x array to give errorbar plots the correct x coordinate
  672. for g = 1:NumberOfFreqClassificationGroups
  673. errorbarx(:,g)=1:AnalysisRegions.numRegions;
  674. end
  676. %%{
  677. %error bars = range of the data
  678. maxValy = maxFreqInFreqGroup-meanFreqInFreqGroup; maxVal = max(maxValy(:)); %find the maximum frequency
  679. minValy = meanFreqInFreqGroup-minFreqInFreqGroup; minVal = min(min(minValy(:)),0.9); %find the minimum frequency
  680. %}
  682. %%{
  683. %error bars = standard deviation of the data
  684. maxValy = meanFreqInFreqGroup+sdFreqInFreqGroup; maxVal = max(maxValy(:)); %find the maximum frequency
  685. minValy = meanFreqInFreqGroup-sdFreqInFreqGroup; minVal = min(min(minValy(:)),0.9); %find the minimum frequency
  686. %}
  688. %{
  689. %error bars = standard error of the data
  690. maxValy = meanFreqInFreqGroup+seFreqInFreqGroup; maxVal = max(maxValy(:)); %find the maximum frequency
  691. minValy = meanFreqInFreqGroup-seFreqInFreqGroup; minVal = min(min(minValy(:)),0.9); %find the minimum frequency
  692. %}
  694. %get y-axis limits
  695. if plotProps.classifyPeakTimesAsFrequencies == 1
  696. %Get powers of 10 to set y limits for the frequency plots
  697. vS = num2str(minVal); b = strfind(vS,'0'); %find zeros in the minimum frequency
  698. if length(b)>1
  699. lowestPower = min(1+strfind(vS,'.'), b(2))-strfind(vS,'.'); %get the position after decimal where the minimum value starts (i.e. 1 = 0.1, 2 = 0.01, 3 = 0.001, 4 = 0.0001, etc...)
  700. minVal = 10.^-lowestPower;
  701. else %otherwise, mean - error went negative
  702. lowestPower = 1;
  703. minVal = 10.^-lowestPower;
  704. end
  706. highestPower = ceil(log10(maxVal)); %get the highest power to set upper y limit
  708. else %otherwise classifying by peak times, not frequency
  709. lowestPower= -min(endpowr);
  710. minVal = 0;
  711. highestPower = 0.1+max([stpowr endpowr]);
  712. end
  714. if plotProps.plotPeakLocationClassificationMedian == 1
  715. labelStr = 'Median';
  716. else
  717. labelStr = 'Mean';
  718. end
  720. figure()
  721. %Indicate where injuries/treatments occurred
  722. hold on
  723. if InjuryIndicies.NumberOfInjuriesOrTreatments > 0 %if there was an injury, shade where the injury occurred
  724. for j = 1:InjuryIndicies.NumberOfInjuriesOrTreatments %for all injuries/treatments...
  725. plot([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j)],[10.^(-lowestPower-1) 10.^highestPower],'--k','LineWidth',plotProps.lineWidth/2)
  726. plot([InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[10.^(-lowestPower-1) 10.^highestPower],'--k','LineWidth',plotProps.lineWidth/2)
  727. fill([InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)-InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j) InjuryInds(1,j)+InjuryInds(2,j)],[10.^(-lowestPower-1) 10.^highestPower 10.^highestPower 10.^(-lowestPower-1)],[0.6 0.6 0.6],'facealpha',0.2,'LineStyle','none')
  728. if InjuryInds(3,j) ~= 1 %if injury/treatment is mid-recoding
  729. text(InjuryInds(1,j) ,min(maxVal*1.5,9.5.^highestPower),InjuryIndicies.InjuryLabels{j},'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  730. else %otherwise injury/treatment is before recording start
  731. text(InjuryInds(1,j) ,min(maxVal*1.5,9.5.^highestPower),[InjuryIndicies.InjuryLabels{j} ' Before Recording Start'],'FontSize',plotProps.FigureFontSize*0.8,'HorizontalAlignment', 'center')
  732. end
  733. end
  734. end
  735. %Now plot mean plus/minus standard deviation in the frequencies of each group
  736. %h = errorbar(errorbarx,meanFreqInFreqGroup',minValy',maxValy','LineWidth',floor(plotProps.lineWidth*3/4));
  737. if plotProps.plotPeakLocationClassificationMedian == 1
  738. h = plot(errorbarx,medianFreqInFreqGroup','.-','LineWidth',floor(plotProps.lineWidth*3/4));
  739. else
  740. h = plot(errorbarx,meanFreqInFreqGroup','.-','LineWidth',floor(plotProps.lineWidth*3/4));
  741. end
  742. hold off
  743. %for c = 1:size(h,2)
  744. for c = 1:size(h,1)
  745. h(c).Color = cc(c,:);
  746. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  747. end
  748. set(gca,'XTick',1:AnalysisRegions.numRegions)
  749. if isfield(AnalysisRegions,'RegionLabels')
  750. set(gca,'XTickLabel',AnalysisRegions.RegionLabels)
  751. end
  752. xlim([0.5 AnalysisRegions.numRegions+0.5])
  753. xlabel('Recording Region')
  754. ylim([minVal 10.^highestPower])
  755. if plotProps.classifyPeakTimesAsFrequencies == 1
  756. ylabel([labelStr ' Frequency of Peaks in Each Band (Hz)'])
  757. title('Peak Frequencies Through Time')
  758. else
  759. ylabel([labelStr ' Peak Times in Each Group (s)'])
  760. title('Peak Locations Through Time')
  761. end
  762. legend(h,freqClassificationNames)
  763. set(gca,'FontSize',plotProps.FigureFontSize)
  764. box on
  765. set(gca,'YTick',10.^[-lowestPower:highestPower])
  766. ax = gca;
  767. ax.YGrid = 'on';
  768. set(gca,'YScale','log')
  772. %% Plot Classification Changes
  774. %Plot group switching dynamics in different recording zones
  775. cc = pink(numZones+2); %set plot colors
  777. figure()
  778. for g = 1:NumberOfleaderClassificationGroups %for each group...
  779. q = zeros(numZones,length(ZoneHistograms.LeaderFollower.Zone1));
  780. for zz = 1:numZones %cycle through zones...
  781. zoneName = strcat('Zone',num2str(zz)); %get the zone name to access the data structure
  782. histArray = ZoneHistograms.LeaderFollower.(zoneName); %for the zone, get the probabilities of being in current group, given past group
  784. q(zz,:) = ZoneHistograms.LeaderFollower.(zoneName)(g,:);
  786. end %end cycle through zones
  787. q = q';
  789. subplot(1,NumberOfleaderClassificationGroups,g)
  790. hold on
  791. h=plot(1:NumberOfleaderClassificationGroups, q,'-','MarkerSize',plotProps.markerSize,'LineWidth',plotProps.lineWidth);
  792. hold off
  793. for c = 1:size(h,1)
  794. h(c).Color = cc(c,:);
  795. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  796. end
  797. box on
  798. grid on
  799. title(['Previously Classified as ' leaderClassificationNames{g}])
  800. ylim([0 1])
  801. xlim([0 NumberOfleaderClassificationGroups+1])
  802. set(gca,'YTick',0:0.1:1)
  803. legend(groupOrder);
  804. ylabel('Probability of Current Classification')
  805. set(gca,'XTick',1:NumberOfleaderClassificationGroups,'XTickLabel',leaderClassificationNames)
  806. xlabel('Current Classification')
  807. set(gca,'FontSize',plotProps.FigureFontSize*0.7)
  809. end
  812. %Plot peak count switching dynamics in different recording zones
  813. c = min(10,NumberOfpkCountClassificationGroups);
  814. r = ceil(NumberOfpkCountClassificationGroups/c);
  817. figure()
  818. t=tiledlayout(r,c,'Padding','none','TileSpacing','compact','Padding','compact');
  819. for g = 1:NumberOfpkCountClassificationGroups %for each group...
  820. q = zeros(numZones,length(ZoneHistograms.PeakCount.Zone1));
  821. for zz = 1:numZones %cycle through zones...
  822. zoneName = strcat('Zone',num2str(zz)); %get the zone name to access the data structure
  823. histArray = ZoneHistograms.PeakCount.(zoneName); %for the zone, get the probabilities of being in current group, given past group
  825. q(zz,:) = ZoneHistograms.PeakCount.(zoneName)(g,:);
  827. end %end cycle through zones
  828. q = q';
  830. nexttile
  831. hold on
  832. h=plot(1:NumberOfpkCountClassificationGroups, q,'-','LineWidth',plotProps.lineWidth);
  833. hold off
  834. for c = 1:size(h,1)
  835. h(c).Color = cc(c,:);
  836. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  837. end
  838. box on
  839. grid on
  840. title(['Previous Peak Count = ' pkCountNames{g}])
  841. ylim([0 1])
  842. xlim([0 NumberOfpkCountClassificationGroups+1])
  843. set(gca,'YTick',0:0.1:1)
  844. legend(groupOrder);
  845. set(gca,'XTick',1:NumberOfpkCountClassificationGroups,'XTickLabel',pkCountNames)
  846. set(gca,'FontSize',plotProps.FigureFontSize*0.7)
  848. end
  849. ylabel(t,'Probability of Current Peak Count','FontSize',plotProps.FigureFontSize)
  850. xlabel(t,'Current Peak Count','FontSize',plotProps.FigureFontSize)
  853. %Plot uniformity switching dynamics in different recording zones
  854. c = min(10,NumberOfunifClassificationGroups);
  855. r = ceil(NumberOfunifClassificationGroups/c);
  857. figure()
  858. t=tiledlayout(r,c,'Padding','none','TileSpacing','compact','Padding','compact');
  859. for g = 1:NumberOfunifClassificationGroups %for each group...
  860. q = zeros(numZones,length(ZoneHistograms.Unif.Zone1));
  861. for zz = 1:numZones %cycle through zones...
  862. zoneName = strcat('Zone',num2str(zz)); %get the zone name to access the data structure
  863. histArray = ZoneHistograms.Unif.(zoneName); %for the zone, get the probabilities of being in current group, given past group
  865. q(zz,:) = ZoneHistograms.Unif.(zoneName)(g,:);
  867. end %end cycle through zones
  868. q = q';
  870. nexttile
  871. hold on
  872. h=plot(1:NumberOfunifClassificationGroups, q,'-','LineWidth',plotProps.lineWidth);
  873. hold off
  874. for c = 1:size(h,1)
  875. h(c).Color = cc(c,:);
  876. h(c).Marker = markerNmArr{c}; h(c).MarkerSize = markerSzArr(c);
  877. end
  878. box on
  879. grid on
  880. title(['Previously ' unifNames{g}])
  881. ylim([0 1])
  882. xlim([0 NumberOfunifClassificationGroups+1])
  883. set(gca,'YTick',0:0.1:1)
  884. legend(groupOrder);
  885. set(gca,'XTick',1:NumberOfunifClassificationGroups,'XTickLabel',unifNames)
  886. set(gca,'FontSize',plotProps.FigureFontSize*0.7)
  888. end
  889. ylabel(t,'Probability of Current Uniformity','FontSize',plotProps.FigureFontSize)
  890. xlabel(t,'Current Uniformity','FontSize',plotProps.FigureFontSize)
  893. %% Plot a summary of all possible correlogram classifications, and the percent of correlograms that fall into each category
  895. totalCorrelogramCount = size(LFGroupAssignments,1)-sum(isnan(LFGroupAssignments)); %get total number of correlograms
  896. UnifMat = zeros(NumberOfpkCountClassificationGroups,NumberOfleaderClassificationGroups,numZones); %initialize array to store classification combos of uniform correlograms
  897. NonUnifMat = zeros(NumberOfpkCountClassificationGroups,NumberOfleaderClassificationGroups,numZones); %initialize array to store classification combos of nonuniform correlograms
  900. for rr = 1:numZones %for each recording region or zone...
  901. rrInds = Zonez == rr;
  902. for u = 1:NumberOfunifClassificationGroups %for each uniformity classification...
  904. nm = strcat('G',num2str(u)); %get the uniformity group name
  905. UIdx=FindGroups.Unif.(nm)(:,rrInds); %find correlograms classified with the given value (nonuniform or uniform)
  907. for p = 1:NumberOfpkCountClassificationGroups %cycle through peak counts
  908. nm2 = strcat('G',num2str(p)); %get the peak count group name
  909. PIdx = FindGroups.PeakCount.(nm2)(:,rrInds); %find correlograms with the given peak count
  911. for l = 1:NumberOfleaderClassificationGroups %cycle through leader/follower groupings
  912. nm3 = strcat('G',num2str(l)); %get the leader/follower group name
  913. LIdx = FindGroups.LeaderFollower.(nm3)(:,rrInds); %find correlograms with the given leader/follower classification
  915. if u == 1 %check in which uniformity matrix to save the entry
  916. NonUnifMat(p,l,rr) = sum(sum(UIdx&PIdx&LIdx))./sum(totalCorrelogramCount(rrInds));
  917. else %otherwise uniform matrix
  918. UnifMat(p,l,rr) = sum(sum(UIdx&PIdx&LIdx))./sum(totalCorrelogramCount(rrInds));
  919. end %end save the fraction of correlograms with a given classification combo in the matrix
  921. end %end cycle through leader/follower groupings
  923. end %end cycle through peak counts
  925. end %end cycle through uniformity classifications
  927. end %end cycle through zones
  929. %Get x and y coordinates for the plot (these are peak counts and leader/follower classifications, uniformity is plotted as different markers)
  930. x = pkCountGroupBoundaries;
  931. y = 1:NumberOfleaderClassificationGroups;
  933. colormapLength = 256;
  934. cc = [1 1 1; jet(colormapLength)]; %get colormap to shade correlogram percentages
  935. xposshift = 0;%0.1; %set a shift factor for the x position so uniform and uniform points are not plotted on top of eachother
  936. yposshift = 0.07; %set a shift factor for the y position so uniform and uniform points are not plotted on top of eachother
  938. for rr = 1:numZones %for each recording zone
  939. figure()
  940. hold on
  941. for xx = 1:length(x) %for each peak count...
  942. for yy = 1:length(y) %for each leader/follower class...
  943. %Plot uniform correlograms
  944. cmapIdx2 = max(1,1+ceil(colormapLength.*UnifMat(xx,yy,rr))); %shade the point according to % of correlograms with the given classification combo
  945. plot(x(xx)+xposshift,y(yy)+yposshift,'^','Color',cc(cmapIdx2,:),'MarkerSize',20,'MarkerFaceColor',cc(cmapIdx2,:),'MarkerEdgeColor',[0 0 0])
  946. text(x(xx)+xposshift+0.1,y(yy)+yposshift+0.1,[num2str(100.*UnifMat(xx,yy,rr),'%.1f') '%']) %display the % of correlograms with the given classification combo
  947. %Plot nonuniform correlograms
  948. cmapIdx = max(1,1+ceil(colormapLength.*NonUnifMat(xx,yy,rr))); %shade according to % of correlograms with the given classification combo
  949. plot(x(xx)-xposshift,y(yy)-yposshift,'v','MarkerSize',20,'MarkerFaceColor',cc(cmapIdx,:),'MarkerEdgeColor',[0 0 0])
  950. text(x(xx)-xposshift+0.1,y(yy)-yposshift-0.1,[num2str(100.*NonUnifMat(xx,yy,rr),'%.1f') '%']) %display the % of correlograms with the given classification combo
  951. end %end cycle through leader/follower classes
  952. end %end cycle through peak counts
  954. %Plot empty points for the legend
  955. h1=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength)),:),'MarkerEdgeColor',[0 0 0]);
  956. h2=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.9)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.9)),:),'MarkerEdgeColor',[0 0 0]);
  957. h3=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.8)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.8)),:),'MarkerEdgeColor',[0 0 0]);
  958. h4=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.7)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.7)),:),'MarkerEdgeColor',[0 0 0]);
  959. h5=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.6)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.6)),:),'MarkerEdgeColor',[0 0 0]);
  960. h6=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.5)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.5)),:),'MarkerEdgeColor',[0 0 0]);
  961. h7=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.4)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.4)),:),'MarkerEdgeColor',[0 0 0]);
  962. h8=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.3)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.3)),:),'MarkerEdgeColor',[0 0 0]);
  963. h9=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.2)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.2)),:),'MarkerEdgeColor',[0 0 0]);
  964. h10=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0.1)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0.1)),:),'MarkerEdgeColor',[0 0 0]);
  965. h11=plot(NaN,NaN,'s','Color',cc(max(1,round(colormapLength*0)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0)),:),'MarkerEdgeColor',[0 0 0]);
  967. hnu = plot(NaN,NaN,'v','Color',cc(max(1,round(colormapLength*0)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0)),:),'MarkerEdgeColor',[0 0 0]);
  968. hu = plot(NaN,NaN,'^','Color',cc(max(1,round(colormapLength*0)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0)),:),'MarkerEdgeColor',[0 0 0]);
  969. hblank = plot(NaN,NaN,'^','Color',cc(max(1,round(colormapLength*0)),:),'MarkerSize',20,'MarkerFaceColor',cc(max(1,round(colormapLength*0)),:),'MarkerEdgeColor',[1 1 1]);
  971. hold off
  972. %legend([hu,hnu,hblank, h1, h2, h3, h4, h5, h6, h7, h8, h9, h10, h11], {'Uniform Correlograms';'Nonuniform Correlograms';'Color = Percent of all Correlograms:'; '100%';'90%';'80%';'70%';'60%';'50%';'40%';'30%';'20%';'10%';'0%'},'Location','EastOutside')
  973. legend([hu,hnu,h1, h2, h3, h4, h5, h6, h7, h8, h9, h10, h11], {'Uniform Correlograms';'Nonuniform Correlograms';'100% Percent of all Correlograms';'90% Percent of all Correlograms';'80% Percent of all Correlograms';'70% Percent of all Correlograms';'60% Percent of all Correlograms';'50% Percent of all Correlograms';'40% Percent of all Correlograms';'30% Percent of all Correlograms';'20% Percent of all Correlograms';'10% Percent of all Correlograms';'0% Percent of all Correlograms'},'Location','EastOutside')
  974. box on
  975. grid on
  976. xlim([min(x)-10*max(xposshift,yposshift) max(x)+10*max(xposshift,yposshift)])
  977. ylim([min(y)-10*max(xposshift,yposshift) max(y)+10*max(xposshift,yposshift)])
  978. set(gca,'XTick',x,'XTickLabel',pkCountNames)
  979. xlabel('Number of Correlogram Peaks')
  980. set(gca,'YTick',y,'YTickLabel',leaderClassificationNames)
  981. ylabel('Leader/Follower Strength Classification')
  982. title(groupOrder{rr})
  983. set(gca,'FontSize',plotProps.FigureFontSize)
  984. end %end cycle through recording zones
  987. end

summarizeCorrelograms.m at commit 1ce9141, under MIT · at the source

Overview

Authors: Casey Erin Adam1,2, Shatha J Mufti1,2,3, Jhon Martinez1,2,3, Edmond A Rogers1,4, Martina Dalolio1,3, Nikita Krishnan1,2,3, Timothy Beauclair1,2,3, Riyi Shi1,2,3
ORCID iDs: Casey Erin Adam, Shatha J Mufti, Jhon Martinez, Edmond A Rogers, Martina Dalolio, Nikita Krishnan, Timothy Beauclair, Riyi Shi
  1. Center for Paralysis Research, Purdue University, West Lafayette, IN 47907 USA
  2. Department of Basic Medical Sciences, School of Veterinary Medicine, Purdue University, West Lafayette, IN 47907 USA
  3. Weldon School of Biomedical Engineering, Purdue University, West Lafayette, IN 47907 USA
  4. Indiana University School of Medicine, 340 West 10 StreetSuite 6200, Fairbanks HallIndianapolis, IN 46202-3082 USA
Institutions: Purdue University West Lafayette (United States); Indiana University School of Medicine
Journal: Neuroinformatics, volume 24, issue 2, article 24
Dates: received 3 September 2025; accepted 26 January 2026; published online 27 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1007/s12021-026-09770-9 · PMID 42036490 · PMCID PMC13111530 · OpenAlex W7155741402
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: traumatic brain injury (population), methods / tools (subfield)
Methods: Preprocessing, Connectivity, Statistics, Single-unit activity, calcium imaging
Keywords: Raster, Correlogram, Bicuculline, Alkalosis, Injury, Software
MeSH: Action Potentials*, Brain Injuries, Traumatic*, Neurons*, Algorithms, Animals, Bicuculline, Microelectrodes, Signal Processing, Computer-Assisted (* major topic)
Topic: Neuroscience and Neural Engineering (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 89 references in the paper

Abstract

Changes in neuronal network dynamics in response to different treatments or conditions underly brain function and pathology. A multitude of tools, including electroencephalography, voltage or calcium imaging, and microelectrode array (MEA) recordings, exist to record signals from neuronal firing at different length scales. Correlograms are a standard analysis tool to study the relationship between a pair of recorded neuronal signals. Correlogram shape can provide information about signal independence, firing pattern, and firing order. However, such analysis is performed manually and qualitatively, limiting the amount of information gained. To overcome this limitation, a MATLAB algorithm was developed to automate correlogram shape quantification by calculating correlogram uniformity, peak count and location, and area left of zero, which respectively quantify signal independence/dependence, firing pattern, and firing order. Algorithm outputs were validated using three different MEA recordings during which cells were exposed to bicuculline methiodide, pH shock, or impact injury. Algorithm outputs described signaling changes in all three recordings, bridged changes in individual signal pairings to changes in the entire signal population, and agreed with literature studies. Therefore, this algorithm serves as a useful means of automatically quantifying changes in signal dependence, firing pattern, and firing order across time within a single recording, and across different recordings. These features are also common to all recording techniques, and therefore can be compared across different types of recordings. Therefore, the MATLAB algorithm described in this article can help provide insight into how neural network dynamics are altered by different drugs or conditions.

Supplementary Information: The online version contains supplementary material available at 10.1007/s12021-026-09770-9.

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

Repository

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

ceadm/Correlogram-Uniformty-Peak-Count-Area-Left-of-Zero

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 1ce9141968e05c8cf5b48fd1a04ce067c1493c59, 10 June 2026
Languages: MATLAB (18), C/C++ (1), C (1)
Size: 24 files, 20 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (2 files), Signal Processing Toolbox (1 file), Violinplot-Matlab (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
22 files
At the source: github.com/ceadm/Correlogram-Uniformty-Peak-Count-Area-Left-of-Zero.git

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

Data Availability

The MATLAB script is available on GitHub (https://github.com/ceadm/Correlogram-Uniformty-Peak-Count-Area-Left-of-Zero.git) as well as in the supplementary material of this work. Recording files are also provided with the code in the supplementary material of this work.

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 6 keywords, 8 MeSH terms, 80 references.

Cite

This paper

Adam, C. E., Mufti, S. J., Martinez, J., Rogers, E. A., Dalolio, M., Krishnan, N., Beauclair, T., & Shi, R. (2026). Quantifying Correlogram Shape to Analyze Neuronal Firing Dynamics Recorded in TBI-on-a-Chip. Neuroinformatics, 24(2), 24. https://doi.org/10.1007/s12021-026-09770-9

BibTeX

@article{adam2026quantifying,
author = {Adam, Casey Erin and Mufti, Shatha J and Martinez, Jhon and Rogers, Edmond A and Dalolio, Martina and Krishnan, Nikita and Beauclair, Timothy and Shi, Riyi},
title = {{Quantifying Correlogram Shape to Analyze Neuronal Firing Dynamics Recorded in TBI-on-a-Chip}},
journal = {Neuroinformatics},
year = {2026},
month = apr,
volume = {24},
number = {2},
pages = {24},
publisher = {Springer Science+Business Media},
issn = {1539-2791},
doi = {10.1007/s12021-026-09770-9},
url = {https://doi.org/10.1007/s12021-026-09770-9},
pmid = {42036490},
pmcid = {PMC13111530}
}

RIS

TY - JOUR
AU - Adam, Casey Erin
AU - Mufti, Shatha J
AU - Martinez, Jhon
AU - Rogers, Edmond A
AU - Dalolio, Martina
AU - Krishnan, Nikita
AU - Beauclair, Timothy
AU - Shi, Riyi
TI - Quantifying Correlogram Shape to Analyze Neuronal Firing Dynamics Recorded in TBI-on-a-Chip
T2 - Neuroinformatics
J2 - Neuroinformatics
PY - 2026
DA - 2026/04/27
VL - 24
IS - 2
SP - 24
SN - 1539-2791
PB - Springer Science+Business Media
DO - 10.1007/s12021-026-09770-9
UR - https://doi.org/10.1007/s12021-026-09770-9
LA - en
ER -

CSL-JSON

{
"id": "10.1007/s12021-026-09770-9",
"type": "article-journal",
"title": "Quantifying Correlogram Shape to Analyze Neuronal Firing Dynamics Recorded in TBI-on-a-Chip",
"container-title": "Neuroinformatics",
"author": [
{
"family": "Adam",
"given": "Casey Erin"
},
{
"family": "Mufti",
"given": "Shatha J"
},
{
"family": "Martinez",
"given": "Jhon"
},
{
"family": "Rogers",
"given": "Edmond A"
},
{
"family": "Dalolio",
"given": "Martina"
},
{
"family": "Krishnan",
"given": "Nikita"
},
{
"family": "Beauclair",
"given": "Timothy"
},
{
"family": "Shi",
"given": "Riyi"
}
],
"container-title-short": "Neuroinformatics",
"volume": "24",
"issue": "2",
"page": "24",
"DOI": "10.1007/s12021-026-09770-9",
"PMID": "42036490",
"PMCID": "PMC13111530",
"ISSN": "1539-2791",
"publisher": "Springer Science+Business Media",
"URL": "https://doi.org/10.1007/s12021-026-09770-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
27
]
]
}
}

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.1080/26941899.2026.2619222
Neurodatascience: Past, Present, and Future.
Journal: Data science in science
In common: 4 references
[2] doi:10.1093/bioinformatics/btag570 [code]
CASCADE: criticality avalanche spike cross-platform analysis detection engine, a multi-manufacturer MEA bash analysis pipeline.
Journal: Bioinformatics (Oxford, England)
In common: methods / tools, 3 references
[3] doi:10.1371/journal.pcbi.1014615 [code]
Toward reliable machine learning models for neural circuit inference: A diagnostic study of CNNs on spike trains.
Journal: PLoS computational biology
In common: 3 references
[4] doi:10.1016/j.bpr.2026.100278
Information flow in cultured neuronal networks by transfer entropy.
Journal: Biophysical reports
In common: 3 references
[5] doi:10.1002/advs.77857 [code]
Brain Network Dynamics of Local and Global Predictive Processing in Aging.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: Violinplot-Matlab, Signal Processing Toolbox, Statistics and Machine Learning Toolbox
[6] doi:10.1038/s41467-026-76581-6 [code]
Thalamocortical bursts encode reward contingencies and drive associative learning.
Journal: Nature communications
In common: Violinplot-Matlab, Signal Processing Toolbox, Statistics and Machine Learning Toolbox
[7] doi:10.1002/glia.70181 [code]
Female Mice Show Stronger Time-of-Day Modulation of Astrocytic Ca&lt;sup&gt;2+&lt;/sup&gt; Activity in the Sleep-Regulatory Ventrolateral Preoptic Nucleus.
Journal: Glia
In common: Violinplot-Matlab, Signal Processing Toolbox, Statistics and Machine Learning Toolbox
[8] doi:10.1038/s41467-026-75492-w [code]
Place and behavioral modulation of hippocampal neurons during immobility.
Journal: Nature communications
In common: Violinplot-Matlab, Signal Processing Toolbox, Statistics and Machine Learning Toolbox
[9] doi:10.1162/netn.a.554 [code]
The turbulent brain: Modeling vortex interactions for understanding human cognition.
Journal: Network neuroscience (Cambridge, Mass.)
In common: Violinplot-Matlab, Signal Processing Toolbox, Statistics and Machine Learning Toolbox
[10] doi:10.1038/s41467-026-75490-y [code]
Topographically organized dorsal raphe activity modulates forebrain sensory-motor representations and contributes to defensive behaviors.
Journal: Nature communications
In common: Violinplot-Matlab, Signal Processing Toolbox, Statistics and Machine Learning Toolbox

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.