CCG - Compute multiple cross- and auto-correlograms, or cross-covariances
USAGE
[ccg,t,tau,c] = CCG(times,id,<options>)
times times of all events (see NOTE below)
id ID for each event (e.g. unit ID) from 1 to n
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'binSize' bin size in s (default = 0.01)
'duration' duration in s of each xcorrelogram (default = 2)
'nBins' number of bins (default = duration/binSize)
'smooth' smoothing size in bins (0 = no smoothing, default)
'groups' group number (1 or 2) for each event, used to restrict
cross-correlograms to pairs across two groups of events
(see EXAMPLE #2 below)
'mode' 'ccg' or 'ccv' (default = 'ccg')
'alpha' significance level to determine correlated pairs
'totalTime' total recording duration in s (if different from the
default = max(times) - min(times))
=========================================================================
NOTES
The size of the cross-correlograms can be supplied either as a duration
(property 'duration') or as an number of bins (property 'nBins').
OUTPUT
ccg value of cross-correlograms or cross-covariances
dimensions are (nbins,m,n) where m is the number of
reference time series (e.g. reference units) and n the
number of referenced time series (in general m = n,
except when using option 'groups')
t time bins
tau lag times for a each pair (mode 'ccv' only)
c maximum cross-covariance for a each pair (mode 'ccv' only)
NOTE
Parameters 'times', 'id' and 'group' can be obtained using <a href="matlab:help CCGParameters">CCGParameters</a>.
As a special case, when computing the correlograms of spike trains, one
can use the output of <a href="matlab:help GetSpikes">GetSpikes</a> either directly or in combination with
<a href="matlab:help CCGParameters">CCGParameters</a>. See EXAMPLES below.
EXAMPLES
% Auto- and cross-correlograms between all neurons
spikes = GetSpikes('output','numbered');
[ccg,t] = CCG(spikes(:,1),spikes(:,2));
% Only tetrode #1 vs tetrode #2 (e.g. mPFC vs HPC neurons)
pfc = GetSpikes([1 -1],'output','numbered');
hpc = GetSpikes([2 -1],'output','numbered');
[s,ids,groups] = CCGParameters(pfc,hpc,2);
[ccg,t] = CCG(s,ids,'groups',groups);
% Between stimulations and MUA spikes
spikes = GetSpikes;
stimulatios = GetEvents('Stimulation');
d = [spikes(:,1) ones(size(spikes,1)) ; stimulations 2*ones(size(stimulations,1))];
d = sortrows(d);
[ccg,t] = CCG(d(:,1),d(:,2));
% To compute cross-covariances
[ccv,t,tau,C] = CCG(times,ids,'mode','ccv');
SEE
See also CCGParameters, ShortTimeCCG.0001 function [ccg,t,tau,c] = CCG(times,id,varargin) 0002 0003 %CCG - Compute multiple cross- and auto-correlograms, or cross-covariances 0004 % 0005 % USAGE 0006 % 0007 % [ccg,t,tau,c] = CCG(times,id,<options>) 0008 % 0009 % times times of all events (see NOTE below) 0010 % id ID for each event (e.g. unit ID) from 1 to n 0011 % <options> optional list of property-value pairs (see table below) 0012 % 0013 % ========================================================================= 0014 % Properties Values 0015 % ------------------------------------------------------------------------- 0016 % 'binSize' bin size in s (default = 0.01) 0017 % 'duration' duration in s of each xcorrelogram (default = 2) 0018 % 'nBins' number of bins (default = duration/binSize) 0019 % 'smooth' smoothing size in bins (0 = no smoothing, default) 0020 % 'groups' group number (1 or 2) for each event, used to restrict 0021 % cross-correlograms to pairs across two groups of events 0022 % (see EXAMPLE #2 below) 0023 % 'mode' 'ccg' or 'ccv' (default = 'ccg') 0024 % 'alpha' significance level to determine correlated pairs 0025 % 'totalTime' total recording duration in s (if different from the 0026 % default = max(times) - min(times)) 0027 % ========================================================================= 0028 % 0029 % NOTES 0030 % 0031 % The size of the cross-correlograms can be supplied either as a duration 0032 % (property 'duration') or as an number of bins (property 'nBins'). 0033 % 0034 % OUTPUT 0035 % ccg value of cross-correlograms or cross-covariances 0036 % dimensions are (nbins,m,n) where m is the number of 0037 % reference time series (e.g. reference units) and n the 0038 % number of referenced time series (in general m = n, 0039 % except when using option 'groups') 0040 % t time bins 0041 % tau lag times for a each pair (mode 'ccv' only) 0042 % c maximum cross-covariance for a each pair (mode 'ccv' only) 0043 % 0044 % 0045 % NOTE 0046 % 0047 % Parameters 'times', 'id' and 'group' can be obtained using <a href="matlab:help CCGParameters">CCGParameters</a>. 0048 % As a special case, when computing the correlograms of spike trains, one 0049 % can use the output of <a href="matlab:help GetSpikes">GetSpikes</a> either directly or in combination with 0050 % <a href="matlab:help CCGParameters">CCGParameters</a>. See EXAMPLES below. 0051 % 0052 % EXAMPLES 0053 % 0054 % % Auto- and cross-correlograms between all neurons 0055 % spikes = GetSpikes('output','numbered'); 0056 % [ccg,t] = CCG(spikes(:,1),spikes(:,2)); 0057 % 0058 % % Only tetrode #1 vs tetrode #2 (e.g. mPFC vs HPC neurons) 0059 % pfc = GetSpikes([1 -1],'output','numbered'); 0060 % hpc = GetSpikes([2 -1],'output','numbered'); 0061 % [s,ids,groups] = CCGParameters(pfc,hpc,2); 0062 % [ccg,t] = CCG(s,ids,'groups',groups); 0063 % 0064 % % Between stimulations and MUA spikes 0065 % spikes = GetSpikes; 0066 % stimulatios = GetEvents('Stimulation'); 0067 % d = [spikes(:,1) ones(size(spikes,1)) ; stimulations 2*ones(size(stimulations,1))]; 0068 % d = sortrows(d); 0069 % [ccg,t] = CCG(d(:,1),d(:,2)); 0070 % 0071 % % To compute cross-covariances 0072 % [ccv,t,tau,C] = CCG(times,ids,'mode','ccv'); 0073 % 0074 % SEE 0075 % 0076 % See also CCGParameters, ShortTimeCCG. 0077 0078 % Copyright (C) 2012-2013 by Michaël Zugaro, Marie Goutierre 0079 % 0080 % This program is free software; you can redistribute it and/or modify 0081 % it under the terms of the GNU General Public License as published by 0082 % the Free Software Foundation; either version 3 of the License, or 0083 % (at your option) any later version. 0084 0085 0086 % Default values 0087 d = 2; 0088 duration = []; 0089 binSize = 0.01; 0090 nBins = []; 0091 smooth = 0; 0092 groups = []; 0093 mode = 'ccg'; 0094 alpha = 0.05; 0095 0096 % Check parameters 0097 if nargin < 2, 0098 error('Incorrect number of parameters (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0099 end 0100 if ~isdvector(times), 0101 error('Parameter ''times'' is not a real-valued vector (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0102 end 0103 if ~isdscalar(id) && ~isdvector(id), 0104 error('Parameter ''id'' is not a real-valued scalar or vector (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0105 end 0106 if ~isdscalar(id) && length(times) ~= length(id), 0107 error('Parameters ''times'' and ''id'' have different lengths (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0108 end 0109 id = id(:); 0110 times = times(:); 0111 totalTime = max(times)-min(times); 0112 0113 % Parse parameter list 0114 for i = 1:2:length(varargin), 0115 if ~ischar(varargin{i}), 0116 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help CCG">CCG</a>'' for details).']); 0117 end 0118 switch(lower(varargin{i})), 0119 case 'binsize', 0120 binSize = varargin{i+1}; 0121 if ~isdscalar(binSize,'>0'), 0122 error('Incorrect value for property ''binSize'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0123 end 0124 case 'nbins', 0125 nBins = varargin{i+1}; 0126 if ~isiscalar(nBins,'>0'), 0127 error('Incorrect value for property ''nBins'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0128 end 0129 case 'duration', 0130 duration = varargin{i+1}; 0131 if ~isdscalar(duration,'>0'), 0132 error('Incorrect value for property ''duration'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0133 end 0134 case 'smooth', 0135 smooth = varargin{i+1}; 0136 if ~isdscalar(smooth,'>=0'), 0137 error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0138 end 0139 case 'groups', 0140 groups = varargin{i+1}; 0141 if ~isempty(groups) && ~isdvector(groups) && length(times) ~= length(groups) 0142 error('Incorrect value for property ''groups'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0143 end 0144 case 'alpha', 0145 alpha = varargin{i+1}; 0146 if ~isdscalar(alpha,'>0'), 0147 error('Incorrect value for property ''alpha'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0148 end 0149 case 'mode', 0150 mode = varargin{i+1}; 0151 if ~isastring(mode,'ccg','ccv'), 0152 error('Incorrect value for property ''mode'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0153 end 0154 case 'totaltime', 0155 totalTime = varargin{i+1}; 0156 if ~isdscalar(totalTime,'>0'), 0157 error('Incorrect value for property ''totaltime'' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0158 end 0159 otherwise, 0160 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help CCG">CCG</a>'' for details).']); 0161 end 0162 end 0163 0164 % Determine binSize/duration 0165 if isempty(nBins), 0166 if isempty(duration), 0167 duration = d; 0168 end 0169 else 0170 if isempty(duration), 0171 duration = nBins*binSize; 0172 elseif duration ~= binSize*nBins, 0173 error('Incompatible ''duration'' and ''nBins'' parameters (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0174 end 0175 end 0176 0177 tau = []; 0178 c = []; 0179 0180 % Number of IDs, number of bins, etc. 0181 if length(id) == 1, 0182 id = ones(length(times),1); 0183 nIDs = 1; 0184 else 0185 nIDs = length(unique(id)); 0186 end 0187 if nIDs ~= max(id), 0188 error('Incorrect IDs (type ''help <a href="matlab:help CCG">CCG</a>'' for details).'); 0189 end 0190 halfBins = round(duration/binSize/2); 0191 nBins = 2*halfBins+1; 0192 t = (-halfBins:halfBins)'*binSize; 0193 0194 if length(times) <= 1, 0195 return 0196 end 0197 0198 % Sort events in time and compute CCGs 0199 [times,i] = sort(times); 0200 id = id(i); 0201 if ~isempty(groups), 0202 groups = groups(i); 0203 end 0204 counts = CCGEngine(times,id,binSize,halfBins); 0205 0206 % Reshape the results 0207 n = max(id); 0208 counts = reshape(counts,[nBins n n]); 0209 if n < nIDs, 0210 counts(nBins,nIDs,nIDs) = 0; 0211 end 0212 0213 % Restrict the results to inter-group CCGs if requested 0214 if ~isempty(groups), 0215 group1 = unique(id(groups == 1)); 0216 group2 = unique(id(groups == 2)); 0217 nGroup1 = length(group1); 0218 nGroup2 = length(group2); 0219 ccg = zeros(nBins,nGroup1,nGroup2); 0220 for i = 1:nGroup1, 0221 for j = 1:nGroup2, 0222 ccg(:,i,j) = Smooth(flipud(counts(:,group1(i),group2(j))),smooth); 0223 end 0224 end 0225 else 0226 ccg = zeros(nBins,nIDs,nIDs); 0227 % Compute corr(A,B) for each unique unordered pair (A,B) 0228 for g1 = 1:nIDs, 0229 for g2 = g1:nIDs, 0230 ccg(:,g1,g2) = Smooth(flipud(counts(:,g1,g2)),smooth); 0231 end 0232 end 0233 % corr(B,A) and corr(B,A) symmetric 0234 for g1 = 1:nIDs, 0235 for g2 = 1:g1-1, 0236 ccg(:,g1,g2) = flipud(squeeze(ccg(:,g2,g1))); 0237 end 0238 end 0239 end 0240 0241 0242 if strcmp(mode,'ccv'), 0243 0244 % Determine mean event rate for each ID 0245 eventRate = zeros(nIDs,1); 0246 for i = 1:nIDs, 0247 eventRate(i) = sum(id==i)/totalTime; 0248 end 0249 0250 % Determine standardized cross-covariances 0251 ccv = zeros(size(ccg)); 0252 tau = zeros(size(ccg,2),size(ccg,3)); 0253 c = zeros(size(ccg,2),size(ccg,3)); 0254 0255 nPairs = size(ccg,2)*size(ccg,3); 0256 disp(['# pairs: ' int2str(nPairs)]); 0257 0258 threshold = sqrt(2)*erfinv(1-(alpha/length(t))); 0259 0260 for i = 1:size(ccg,2), 0261 for j = 1:size(ccg,3), 0262 0263 % Compute and normalize CCVs from CCGs 0264 if ~isempty(groups), 0265 rate = eventRate(group1(i))*eventRate(group2(j)); 0266 else 0267 rate = eventRate(i)*eventRate(j); 0268 end 0269 ccv(:,i,j) = sqrt((binSize*totalTime)/rate) * (ccg(:,i,j)/(binSize*totalTime)-rate); 0270 0271 % Smooth with a 3-bin boxcar 0272 data = ccv(:,i,j); 0273 top = flipud(data(1:size(ccg,1),:)); 0274 bottom = flipud(data(end-size(ccg,1)+1:end,:)); 0275 data = [top;data;bottom]; 0276 data = filter([1 1 1],3,data); 0277 n = size(data,1); 0278 d = n - size(ccg,1); 0279 start = d/2+1; 0280 stop = start + size(ccg,1) - 1; 0281 ccv(:,i,j) = data(start:stop); 0282 0283 % Find the peak lag time and value 0284 [~,maxIndex] = max(ccv(:,i,j)); 0285 tau(i,j) = median(t(maxIndex)); 0286 c(i,j) = median(ccv(maxIndex,i,j)); 0287 % Previous version of the code (discard?) 0288 % [~,maxIndex] = max(ccv(:,i,j)); 0289 % tau(i,j) = t(maxIndex); 0290 % c(i,j) = median(ccv(max(1,maxIndex-3):min(end,maxIndex+3),i,j)); 0291 0292 % Keep only significantly correlated pairs 0293 if ~any(abs(ccv(:,i,j))>threshold), 0294 tau(i,j) = NaN; 0295 end 0296 0297 end 0298 end 0299 0300 nCorrelatedPairs = sum(~isnan(tau(:))); 0301 disp(['# significantly correlated pairs: ' int2str(nCorrelatedPairs)]); 0302 0303 ccg = ccv; 0304 0305 end