FunctionalClustering - Determine spike train similarity tree using FCA analysis.
Determine which spike trains are most similar in a neuronal population. This
uses the Functional Clustering Algorithm (FCA) of Feldt et al. (2009).
For each pair of neurons i and j, we first measure the distance (time elapsed)
between each spike of i and the closest spike of j, and use the average of these
distances as an estimate of the similarity between i and j. We then swap the
two neurons and measure the similarity between j and i. The average minimal
distance (AMD) is defined as the mean between the two similarity values. Here,
we use a normalizing factor so that a pair of independent Poisson processes
has an AMD of 1.
To define the scaled significance, the spike trains are repeatedly jittered by
a random amount, and the distribution of jittered spike trains is determined.
The AMD is centered and normalized using the median and 95% percentile of this
distribution. Thus, a scaled significance of 1 corresponds to a p value of 0.05.
Finally, the similarity tree is computed recursively by finding the two neurons
with the highest AMD, grouping them into a new cluster, and repeating until all
neurons have been grouped together.
USAGE
[d,s,linkage] = FunctionalClustering(spikes,<options>)
spikes two-column matrix of timestamps and unit IDs,
provided by <a href="matlab:help GetSpikeTimes">GetSpikeTimes</a> using 'output' = 'numbered'
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'jitter' jitter standard deviation, in s (default = 5)
'nJitters' number of jittered spike trains (default = 1000)
=========================================================================
OUTPUT
d normalized average minimal distance for each unit pair
s scaled significance for each unit pair
linkage network connectivity tree
SEE
See also PlotLinkage.0001 function [d,s,linkage] = FunctionalClustering(spikes,varargin) 0002 0003 %FunctionalClustering - Determine spike train similarity tree using FCA analysis. 0004 % 0005 % Determine which spike trains are most similar in a neuronal population. This 0006 % uses the Functional Clustering Algorithm (FCA) of Feldt et al. (2009). 0007 % 0008 % For each pair of neurons i and j, we first measure the distance (time elapsed) 0009 % between each spike of i and the closest spike of j, and use the average of these 0010 % distances as an estimate of the similarity between i and j. We then swap the 0011 % two neurons and measure the similarity between j and i. The average minimal 0012 % distance (AMD) is defined as the mean between the two similarity values. Here, 0013 % we use a normalizing factor so that a pair of independent Poisson processes 0014 % has an AMD of 1. 0015 % 0016 % To define the scaled significance, the spike trains are repeatedly jittered by 0017 % a random amount, and the distribution of jittered spike trains is determined. 0018 % The AMD is centered and normalized using the median and 95% percentile of this 0019 % distribution. Thus, a scaled significance of 1 corresponds to a p value of 0.05. 0020 % 0021 % Finally, the similarity tree is computed recursively by finding the two neurons 0022 % with the highest AMD, grouping them into a new cluster, and repeating until all 0023 % neurons have been grouped together. 0024 % 0025 % USAGE 0026 % 0027 % [d,s,linkage] = FunctionalClustering(spikes,<options>) 0028 % 0029 % spikes two-column matrix of timestamps and unit IDs, 0030 % provided by <a href="matlab:help GetSpikeTimes">GetSpikeTimes</a> using 'output' = 'numbered' 0031 % <options> optional list of property-value pairs (see table below) 0032 % 0033 % ========================================================================= 0034 % Properties Values 0035 % ------------------------------------------------------------------------- 0036 % 'jitter' jitter standard deviation, in s (default = 5) 0037 % 'nJitters' number of jittered spike trains (default = 1000) 0038 % ========================================================================= 0039 % 0040 % OUTPUT 0041 % 0042 % d normalized average minimal distance for each unit pair 0043 % s scaled significance for each unit pair 0044 % linkage network connectivity tree 0045 % 0046 % SEE 0047 % 0048 % See also PlotLinkage. 0049 0050 % Copyright (C) 2015 by Ralitsa Todorova, Michaël Zugaro 0051 % 0052 % This program is free software; you can redistribute it and/or modify 0053 % it under the terms of the GNU General Public License as published by 0054 % the Free Software Foundation; either version 3 of the License, or 0055 % (at your option) any later version. 0056 0057 % Defaults 0058 jitterSD = 5; 0059 nJitters = 1000; 0060 0061 % Check parameter 0062 if ~isdmatrix(spikes,'@2') || ~isivector(spikes(:,2)), 0063 error('Incorrect spikes (type ''help <a href="matlab:help FunctionalClustering">FunctionalClustering</a>'' for details).'); 0064 end 0065 0066 % Parse parameter list 0067 for i = 1:2:length(varargin), 0068 if ~ischar(varargin{i}), 0069 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help FunctionalClustering">FunctionalClustering</a>'' for details).']); 0070 end 0071 switch(lower(varargin{i})), 0072 case 'jitter', 0073 jitterSD = varargin{i+1}; 0074 if ~isdscalar(jitterSD,'>0'), 0075 error('Incorrect value for property ''jitter'' (type ''help <a href="matlab:help FunctionalClustering">FunctionalClustering</a>'' for details).'); 0076 end 0077 case 'njitters', 0078 nJitters = varargin{i+1}; 0079 if ~isiscalar(nJitters,'>0'), 0080 error('Incorrect value for property ''nJitters'' (type ''help <a href="matlab:help FunctionalClustering">FunctionalClustering</a>'' for details).'); 0081 end 0082 otherwise, 0083 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help FunctionalClustering">FunctionalClustering</a>'' for details).']); 0084 end 0085 end 0086 0087 spikes = sortrows(spikes); 0088 t = spikes(:,1); 0089 id = spikes(:,2); 0090 n = max(id); 0091 T = max(t) - min(t); 0092 % Make sure jittered spikes remain within the time window if this is small 0093 jitterSD = min([jitterSD T/10]); 0094 0095 % The following two matrices have sizes 2n x 2n (rather than n x n) because this will speed up 0096 % computations in the second part of the algorithm, but only the upper n x n square is actually 0097 % used to compute the AMD and scaled significance of the unit pairs 0098 AMD = nan(2*n); 0099 scaledSignificance = nan(2*n); 0100 linkage = nan(1); 0101 jitteredAMD = nan(n,n,nJitters); 0102 0103 % Compute average firing rate for each unit 0104 for i = 1:n, 0105 f(i,1) = sum(id==i)/T; 0106 end 0107 0108 % Compute normalized average minimal distances for each pair of units 0109 for i = 1:n, 0110 for j = i+1:n, 0111 AMD(i,j) = ComputeAMD(spikes,f,i,j); 0112 end 0113 end 0114 d = AMD(1:n,1:n); 0115 0116 % Compute normalized average minimal distances for each pair of jittered units (for each jitter) 0117 for jitter = 1:nJitters, 0118 jittered = sortrows([t+jitterSD*randn(size(t)) id]); 0119 for i = 1:n, 0120 for j = i+1:n, 0121 jitteredAMD(i,j,jitter) = ComputeAMD(jittered,f,i,j); 0122 end 0123 end 0124 end 0125 0126 % Center and reduce to get the scaled significance 0127 medians = quantile(jitteredAMD,0.5,3); 0128 scale = quantile(jitteredAMD,0.05,3) - medians; 0129 s = (AMD(1:n,1:n)-medians)./scale; 0130 0131 % Compute linkage only if required (time consuming) 0132 if nargout < 3, return; end 0133 0134 % Merge most similar spike trains and update matrix 0135 scaledSignificance(1:n,1:n) = s; 0136 linkage = nan(numel(unique(id))-1,4); 0137 iteration = 0; 0138 jitteredAMD = nan(n,nJitters); 0139 while numel(unique(id))>1, 0140 iteration = iteration+1; 0141 % Find max scaled significance 0142 ss = max(scaledSignificance(:)); 0143 [i,j] = find(scaledSignificance==ss,1,'first'); 0144 % Drop the corresponding line and column 0145 scaledSignificance(:,[i j]) = nan; 0146 scaledSignificance([i j],:) = nan; 0147 % Update connectivity tree 0148 linkage(iteration,:) = [i j ss AMD(i,j)]; 0149 % Merge these units into a new cluster 0150 n = n+1; 0151 spikes(id==i|id==j,2) = n; 0152 id = spikes(:,2); 0153 f(n,1) = sum(id==n)/T; 0154 0155 % update AMD and scaled significance 0156 for jitter = 1:nJitters, 0157 jittered = sortrows([t+jitterSD*randn(size(t)) id]); 0158 for i = 1:n-1, 0159 jitteredAMD(i,jitter) = ComputeAMD(jittered,f,i,n); 0160 end 0161 end 0162 for i=1:n-1, 0163 AMD(i,n) = ComputeAMD(spikes,f,i,n);; 0164 medians = quantile(jitteredAMD(i,:),0.5,2); 0165 scale = quantile(jitteredAMD(i,:),0.05,2) - medians; 0166 scaledSignificance(i,n) = (AMD(i,n)-medians)./scale; 0167 end 0168 end 0169 0170 % ------------------------------- Helper function ------------------------------- 0171 0172 function AMD = ComputeAMD(spikes,f,i,j) 0173 0174 AMD = nan; 0175 spikes_ij = spikes(ismember(spikes(:,2),[i j]),:); 0176 if numel(unique(spikes_ij(:,2))) < 2, return; end; 0177 % Distances between successive spikes, whatever unit they belong to 0178 D = [1000 1000; diff(spikes_ij); 1000 1000]; 0179 distance = D(:,1); 0180 different = logical(D(:,2)); % same or different unit? 0181 % Distance to previous spike of the other unit 0182 previous = CumSum(distance,different); 0183 % Distance to next spike of the other unit (flip vectors upside down) 0184 next = flipud(CumSum(flipud(distance),flipud(different))); 0185 % Minimum between the two (for each spike) 0186 m = min([[NaN; previous(2:end-1)] [next(2:end-1); NaN]],[],2); 0187 % (Normalized) average minimal distance (mean of average minimal distances from i to j and from j to i) 0188 d_ij = mean(m(spikes_ij(:,2)==i)); 0189 d_ji = mean(m(spikes_ij(:,2)==j)); 0190 AMD = mean([d_ij d_ji]) / ((f(i)+f(j))/(4*f(i)*f(j)));