ActivityTemplates

ActivityTemplates - Compute activity templates from PCA of spike trains.

 Computes the templates for the component activation analysis described in
 Peyrache et al (2009). These templates can then be tested on different data
 sets using <a href="matlab:help ReactivationStrength">ReactivationStrength</a>. Time bins can be automatically determined
 using a fixed bin size, or provided as an explicit list (e.g. computed using
 theta phases).

  USAGE

    [templates,correlations,eigenvalues] = ActivityTemplates(spikes,<options>)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'bins'        list of [start stop] for all bins
     'binSize'     bin size in s (default = 0.050)
    =========================================================================

  OUTPUT

    templates      3D array of template matrices (dimension 3 is template #,
                   ordered in descending order of corresponding eigenvalue)
    correlations   spike count correlation matrix
    eigenvalues    significant eigenvalues, listed in descending order

  SEE

    See also ReactivationStrength.

0001 function [templates,correlations,eigenvalues,eigenvectors] = ActivityTemplates(spikes,varargin)
0002 
0003 %ActivityTemplates - Compute activity templates from PCA of spike trains.
0004 %
0005 % Computes the templates for the component activation analysis described in
0006 % Peyrache et al (2009). These templates can then be tested on different data
0007 % sets using <a href="matlab:help ReactivationStrength">ReactivationStrength</a>. Time bins can be automatically determined
0008 % using a fixed bin size, or provided as an explicit list (e.g. computed using
0009 % theta phases).
0010 %
0011 %  USAGE
0012 %
0013 %    [templates,correlations,eigenvalues] = ActivityTemplates(spikes,<options>)
0014 %
0015 %    =========================================================================
0016 %     Properties    Values
0017 %    -------------------------------------------------------------------------
0018 %     'bins'        list of [start stop] for all bins
0019 %     'binSize'     bin size in s (default = 0.050)
0020 %    =========================================================================
0021 %
0022 %  OUTPUT
0023 %
0024 %    templates      3D array of template matrices (dimension 3 is template #,
0025 %                   ordered in descending order of corresponding eigenvalue)
0026 %    correlations   spike count correlation matrix
0027 %    eigenvalues    significant eigenvalues, listed in descending order
0028 %
0029 %  SEE
0030 %
0031 %    See also ReactivationStrength.
0032 
0033 % Copyright (C) 2016-2018 by Michaël Zugaro, Ralitsa Todorova
0034 %
0035 % This program is free software; you can redistribute it and/or modify
0036 % it under the terms of the GNU General Public License as published by
0037 % the Free Software Foundation; either version 3 of the License, or
0038 % (at your option) any later version.
0039 
0040 % Defaults
0041 bins = [];
0042 defaultBinSize = 0.050;
0043 binSize = [];
0044 overlap = [];
0045 
0046 % Check number of parameters
0047 if nargin < 1,
0048     error('Incorrect number of parameters (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).');
0049 end
0050 % Check parameter sizes
0051 if ~isdmatrix(spikes,'@2'),
0052     error('Parameter ''spikes'' is not a Nx2 matrix (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).');
0053 end
0054 % Parse parameter list
0055 for i = 1:2:length(varargin),
0056     if ~ischar(varargin{i}),
0057         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).']);
0058     end
0059     switch(lower(varargin{i})),
0060         case 'binsize',
0061             binSize = varargin{i+1};
0062             if ~isdscalar(binSize,'>0'),
0063                 error('Incorrect value for property ''binSize'' (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).');
0064             end
0065         case 'bins',
0066             bins = varargin{i+1};
0067             if ~isdmatrix(bins,'@2'),
0068                 error('Incorrect value for property ''bins'' (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).');
0069             end
0070         otherwise,
0071             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).']);
0072         end
0073 end
0074 
0075 % Options binSize and bins are incompatible
0076 if ~isempty(binSize) && ~isempty(bins),
0077     error('Parameters ''binSize'' and ''bins'' are incompatible (type ''help <a href="matlab:help ActivityTemplates">ActivityTemplates</a>'' for details).');
0078 end
0079 if isempty(binSize) && isempty(bins),
0080     binSize = defaultBinSize;
0081 end
0082 
0083 templates = NaN;
0084 correlations = NaN;
0085 eigenvalues = NaN;
0086 
0087 nUnits = max(spikes(:,2));
0088 if isempty(nUnits), return; end
0089 
0090 %% Bin spikes
0091 spikes = sortrows(spikes,1);
0092 id = spikes(:,2);
0093 
0094 % Shift spike times to start at 0, and list bins unless explicitly provided
0095 if isempty(bins),
0096     spikes(:,1) = spikes(:,1) - spikes(1,1);
0097     bins = (0:binSize:(spikes(end,1)-binSize))';
0098     bins(:,2) = bins+binSize;
0099 else
0100     m = min([min(spikes(:,1)) min(bins(:))]);
0101     spikes(:,1) = spikes(:,1) - m;
0102     bins = bins - m;
0103 end
0104 
0105 % Create spike count matrix
0106 nBins = size(bins,1);
0107 if isempty(nBins), return; end
0108 n = zeros(nBins,nUnits);
0109 for unit = 1:nUnits,
0110     n(:,unit) = CountInIntervals(spikes(id==unit,1),bins);
0111 end
0112 
0113 %% Create correlation matrix
0114 n = zscore(n);
0115 correlations = (1/(nBins-1))*n'*n;
0116 
0117 % Compute eigenvalues/vectors and sort in descending order
0118 [eigenvectors,eigenvalues] = eig(correlations);
0119 [eigenvalues,i] = sort(diag(eigenvalues),'descend');
0120 eigenvectors = eigenvectors(:,i);
0121 
0122 %% Keep only significant eigenvalues and compute templates
0123 
0124 q = nBins/nUnits;
0125 if q < 1,
0126     warning('Not enough time bins to determine significant templates');
0127     eigenvalues = NaN;
0128 end
0129 lambdaMax = (1+sqrt(1/q))^2;
0130 significant = eigenvalues>lambdaMax;
0131 eigenvectors = eigenvectors(:,significant);
0132 templates = zeros(nUnits,nUnits,sum(significant));
0133 for i = 1:sum(significant),
0134     templates(:,:,i) = eigenvectors(:,i)*eigenvectors(:,i)';
0135     templates(:,:,i) = templates(:,:,i) - diag(diag(templates(:,:,i))); % remove the diagonal
0136 end