ReactivationStrength

ReactivationStrength - Assess reactivation strength of cell assemblies.

 Estimates how similar spiking activity is to given templates, across time.
 Time bins can be automatically determined using a fixed bin size, or provided
 as an explicit list (e.g. computed using theta phases).

  USAGE

    strength = ReactivationStrength(spikes,templates,<options>)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'bins'        list of [start stop] for all bins
     'binSize'     bin size in s (default = 0.050)
     'overlap'     overlap between successive bins (default = binSize/2)
     'step'        step between successive bins (default = binSize/2)
    =========================================================================

  OUTPUT

    strength       reactivation strength across time (if time bins are not
                   explicitly provided, the first bin is centered on the
                   first spike)

  SEE

    See also ActivityTemplates.

0001 function strength = ReactivationStrength(spikes,templates,varargin);
0002 
0003 %ReactivationStrength - Assess reactivation strength of cell assemblies.
0004 %
0005 % Estimates how similar spiking activity is to given templates, across time.
0006 % Time bins can be automatically determined using a fixed bin size, or provided
0007 % as an explicit list (e.g. computed using theta phases).
0008 %
0009 %  USAGE
0010 %
0011 %    strength = ReactivationStrength(spikes,templates,<options>)
0012 %
0013 %    =========================================================================
0014 %     Properties    Values
0015 %    -------------------------------------------------------------------------
0016 %     'bins'        list of [start stop] for all bins
0017 %     'binSize'     bin size in s (default = 0.050)
0018 %     'overlap'     overlap between successive bins (default = binSize/2)
0019 %     'step'        step between successive bins (default = binSize/2)
0020 %    =========================================================================
0021 %
0022 %  OUTPUT
0023 %
0024 %    strength       reactivation strength across time (if time bins are not
0025 %                   explicitly provided, the first bin is centered on the
0026 %                   first spike)
0027 %
0028 %  SEE
0029 %
0030 %    See also ActivityTemplates.
0031 
0032 % Copyright (C) 2016-2018 by Michaël Zugaro, Ralitsa Todorova
0033 %
0034 % This program is free software; you can redistribute it and/or modify
0035 % it under the terms of the GNU General Public License as published by
0036 % the Free Software Foundation; either version 3 of the License, or
0037 % (at your option) any later version.
0038 
0039 % Defaults
0040 defaultBinSize = 0.050;
0041 binSize = [];
0042 overlap = [];
0043 step = [];
0044 bins = [];
0045 
0046 for i = 1:2:length(varargin),
0047     if ~ischar(varargin{i}),
0048         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).']);
0049     end
0050     switch(lower(varargin{i})),
0051         case 'binsize',
0052             binSize = varargin{i+1};
0053             if ~isdscalar(binSize,'>0'),
0054                 error('Incorrect value for property ''binSize'' (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).');
0055             end
0056         case 'overlap',
0057             overlap = varargin{i+1};
0058             if ~isdscalar(overlap,'>0'),
0059                 error('Incorrect value for property ''overlap'' (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).');
0060             end
0061         case 'step',
0062             step = varargin{i+1};
0063             if ~isdscalar(step,'>0'),
0064                 error('Incorrect value for property ''step'' (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).');
0065             end
0066         case 'bins',
0067             bins = varargin{i+1};
0068             if ~isdmatrix(bins,'@2'),
0069                 error('Incorrect value for property ''bins'' (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).');
0070             end
0071         otherwise,
0072             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).']);
0073     end
0074 end
0075 
0076 % Option bins is incompatible with binSize, step and overlap
0077 if ~isempty(bins) && (~isempty(binSize) || ~isempty(overlap) || ~isempty(step)) ,
0078     error('Parameter ''bins'' is incompatible with ''binSize'', ''overlap'' and ''step'' (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).');
0079 end
0080 % Unless explicit bins were provided, update binSize and step/overlap
0081 if isempty(bins),
0082     if isempty(binSize),
0083         binSize = defaultBinSize;
0084     end
0085     if isempty(step) && isempty(overlap),
0086         step = binSize/2;
0087     else
0088         if isempty(step),
0089             step = binSize-overlap;
0090         elseif step ~= binSize-overlap,
0091             error('Incompatible ''step'' and ''overlap'' parameters (type ''help <a href="matlab:help ReactivationStrength">ReactivationStrength</a>'' for details).');
0092         end
0093     end
0094 end
0095 
0096 % Shift spikes to start at 0 and bin them
0097 nUnits = size(templates,2);
0098 spikes = spikes(spikes(:,2)<=nUnits,:);
0099 spikes = sortrows(spikes,1);
0100 id = spikes(:,2);
0101 if isempty(bins),
0102     spikes(:,1) = spikes(:,1);
0103     bins = (spikes(1,1):step:(spikes(end,1)-binSize))';
0104     bins(:,2) = bins + binSize;
0105 end
0106 
0107 % Create and compute spike count matrix
0108 nBins = size(bins,1);
0109 dN = zeros(nBins,nUnits);
0110 for unit = 1:nUnits,
0111     dN(:,unit) = CountInIntervals(spikes(id==unit,1),bins);
0112 end
0113 %  dN = squeeze(reshape(dN,1,[],nUnits));
0114 dN = zscore(dN);
0115 
0116 % Compute reactivation strengths
0117 nTemplates = size(templates,3);
0118 strength = zeros(nBins,nTemplates);
0119 for i = 1:nTemplates,
0120     template = templates(:,:,i);
0121     % Set the diagonal to zero to not count coactivation of i and j when i=j
0122     template = template - diag(diag(template));
0123     strength(:,i) = nansum(dN*(template).*dN,2);
0124 end
0125 t = nanmean(bins,2);
0126 strength = [t strength];