A semi-automated MEA spike sorting method for high-throughput assessment of cultured neurons.

1. function [initial_idx_list, HDT_flag] = fix_undersorting(initial_idx_list, spikes, Times, score, refractoryT, params)
2. % Fix undersorting using Hartigan's dip test for bimodality
3. %
4. % INPUTS:
5. % initial_idx_list - Cell array of indices for each cluster after outlier removal
6. % spikes - Spike waveforms matrix (time points × spikes)
7. % Times - Spike timing vector
8. % score - PCA scores for each spike
9. % refractoryT - Refractory period in seconds
10. % params - Parameter structure
11. %
12. % OUTPUTS:
13. % initial_idx_list - Updated cell array of indices after fixing undersorting
14. % HDT_flag - Flag indicating if Hartigan's dip test found bimodality
16. HDT_flag = 0;
17. nboot = 500; % Bootstrap sample size for the dip test
18. count_fix_undersorting = 0;
19. copy_initial_idx_list = initial_idx_list;
20. initial_idx_list = cell(1, length(copy_initial_idx_list)*2); % Pre-allocate larger array
22. for num_unit_initial = 1:length(copy_initial_idx_list)
23. try
24. if length(copy_initial_idx_list{num_unit_initial}) > 10
25. % Get PCA scores for this cluster
26. samplePCA1 = score(copy_initial_idx_list{num_unit_initial}, 1)';
27. samplePCA2 = score(copy_initial_idx_list{num_unit_initial}, 2)';
28. % Perform Hartigan's dip test on the first two PCs
29. [~, p_value1, ~, ~] = HartigansDipSignifTest(samplePCA1, nboot);
30. [~, p_value2, ~, ~] = HartigansDipSignifTest(samplePCA2, nboot);
31. % If either PC shows significant bimodality (p < 0.05)
32. if (p_value1 < 0.05) || (p_value2 < 0.05)
33. try
34. % Perform additional clustering on this unit
35. sample2d = score(copy_initial_idx_list{num_unit_initial}, 1:2);
37. % Use evalclusters to determine optimal K
38. klist = 2:5;
39. myfunc = @(X, K) spectralcluster(X, K);
40. eva = evalclusters(sample2d, myfunc, "DaviesBouldin", 'KList', klist);
41. optimal_K = eva.OptimalK;
43. % Perform spectral clustering with optimal K
44. [cluster_idx, ~] = spectralcluster(sample2d, optimal_K);
45. catch ME
46. if contains(ME.message, 'Invalid data type') || contains(ME.message, 'real array')
47. try
48. % Try k-means as fallback
49. optimal_K = 2; % Default to 2 for fallback
50. [cluster_idx, ~] = kmeans(sample2d, optimal_K, 'Replicates', 2);
51. warning('Used k-means fallback');
52. catch
53. % If even k-means fails, skip this iteration
54. warning('Both spectral clustering and k-means failed');
55. continue;
56. end
57. else
58. rethrow(ME); % Re-throw if it's a different error
59. end
60. end
62. % Split the cluster into K sub-clusters based on optimal K
63. sub_clusters = cell(1, optimal_K);
64. for k = 1:optimal_K
65. sub_clusters{k} = copy_initial_idx_list{num_unit_initial}(cluster_idx == k);
66. end
68. % Validate sub-clusters through pairwise overlap analysis
69. valid_clusters = true(1, optimal_K);
70. merge_groups = cell(1, optimal_K);
71. for k = 1:optimal_K
72. merge_groups{k} = k;
73. end
75. % Check all pairs of sub-clusters for overlap
76. for k1 = 1:optimal_K-1
77. for k2 = k1+1:optimal_K
78. if valid_clusters(k1) && valid_clusters(k2)
79. % Get spike waveforms and times for each sub-cluster
80. spk1 = spikes(:, sub_clusters{k1});
81. t1 = Times(:, sub_clusters{k1});
82. spk2 = spikes(:, sub_clusters{k2});
83. t2 = Times(:, sub_clusters{k2});
85. % Calculate overlap between the two sub-clusters
86. [ovlp1, ovlp2, tdiff] = calculate_overlap(spk1, spk2, t1, t2);
87. ovlp = max(sum(ovlp2 > params.overlap_threshold) / length(ovlp1), sum(ovlp2 > params.overlap_threshold) / length(ovlp2));
89. % Check if the sub-clusters are likely from the same spike
90. sameSpike = (ovlp > params.overlap_threshold) && (tdiff(1) > refractoryT);
92. if sameSpike
93. % Mark for merging - combine indices
94. sub_clusters{k1} = [sub_clusters{k1}; sub_clusters{k2}];
95. valid_clusters(k2) = false;
96. end
97. end
98. end
99. end
101. % Add validated sub-clusters to the output
102. num_valid = sum(valid_clusters);
103. if num_valid > 1
104. % Split detected - add each valid sub-cluster
105. HDT_flag = 1;
106. for k = 1:optimal_K
107. if valid_clusters(k)
108. count_fix_undersorting = count_fix_undersorting + 1;
109. initial_idx_list{count_fix_undersorting} = sub_clusters{k};
110. end
111. end
112. else
113. % All sub-clusters merged back - keep as single unit
114. count_fix_undersorting = count_fix_undersorting + 1;
115. initial_idx_list{count_fix_undersorting} = copy_initial_idx_list{num_unit_initial};
116. end
117. else
118. % Keep as a single unit if no bimodality detected
119. count_fix_undersorting = count_fix_undersorting + 1;
120. initial_idx_list{count_fix_undersorting} = copy_initial_idx_list{num_unit_initial};
121. end
122. end
123. catch
124. % Keep as is if an error occurs
125. fprintf('Error in dip test for unit %d, keeping as is\n', num_unit_initial);
126. count_fix_undersorting = count_fix_undersorting + 1;
127. initial_idx_list{count_fix_undersorting} = copy_initial_idx_list{num_unit_initial};
128. end
129. end
131. % Trim empty cells
132. initial_idx_list = initial_idx_list(~cellfun('isempty', initial_idx_list));
133. end

fix_undersorting.m at commit 9e40131, under MIT · at the source