PhaseMango - Compute phase as a function of spike rate and acceleration.
Compute spike phase as a function of spike rate and acceleration, as in
Harris et al. (2002). The name of the function refers to the aspect of
the resulting phase plot.
USAGE
[map,data] = PhaseMango(spikes,phases,<options>)
spikes spike timestamps
phases phase <a href="matlab:help samples">samples</a> (see <a href="matlab:help Phase">Phase</a>)
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'smooth' smoothing size in bins (0 = no smoothing, default = 2)
'nBins' number of horizontal and vertical bins (default = [50 50])
'limits' [maxRate minAccel maxAccel] (default = [8 -2 2])
'maxGap' time gaps between successive position samples exceeding
this threshold (e.g. undetects) will not be interpolated
(default = 100 ms)
'boundaries' onset and offset for single-lap phase precession can be
determined either automatically based on spike count
('count', default) or using explicit firing field
boundaries ([Xstart Xstop], where each X is in [0..1])
=========================================================================
OUTPUT
map.x x bins
map.y y bins
map.phase phase map
map.count count map
data.rate spike rate for each cycle
data.acceleration spike acceleration for each cycle
data.phase mean spike phase (in radians) for each cycle
SEE
See also Phase, PhasePrecession, PlotPhaseMango.0001 function [map,data] = PhaseMango(spikes,phases,varargin) 0002 0003 %PhaseMango - Compute phase as a function of spike rate and acceleration. 0004 % 0005 % Compute spike phase as a function of spike rate and acceleration, as in 0006 % Harris et al. (2002). The name of the function refers to the aspect of 0007 % the resulting phase plot. 0008 % 0009 % USAGE 0010 % 0011 % [map,data] = PhaseMango(spikes,phases,<options>) 0012 % 0013 % spikes spike timestamps 0014 % phases phase <a href="matlab:help samples">samples</a> (see <a href="matlab:help Phase">Phase</a>) 0015 % <options> optional list of property-value pairs (see table below) 0016 % 0017 % ========================================================================= 0018 % Properties Values 0019 % ------------------------------------------------------------------------- 0020 % 'smooth' smoothing size in bins (0 = no smoothing, default = 2) 0021 % 'nBins' number of horizontal and vertical bins (default = [50 50]) 0022 % 'limits' [maxRate minAccel maxAccel] (default = [8 -2 2]) 0023 % 'maxGap' time gaps between successive position samples exceeding 0024 % this threshold (e.g. undetects) will not be interpolated 0025 % (default = 100 ms) 0026 % 'boundaries' onset and offset for single-lap phase precession can be 0027 % determined either automatically based on spike count 0028 % ('count', default) or using explicit firing field 0029 % boundaries ([Xstart Xstop], where each X is in [0..1]) 0030 % ========================================================================= 0031 % 0032 % OUTPUT 0033 % 0034 % map.x x bins 0035 % map.y y bins 0036 % map.phase phase map 0037 % map.count count map 0038 % 0039 % data.rate spike rate for each cycle 0040 % data.acceleration spike acceleration for each cycle 0041 % data.phase mean spike phase (in radians) for each cycle 0042 % 0043 % SEE 0044 % 0045 % See also Phase, PhasePrecession, PlotPhaseMango. 0046 0047 % Copyright (C) 2004-2012 by Michaël Zugaro 0048 % 0049 % This program is free software; you can redistribute it and/or modify 0050 % it under the terms of the GNU General Public License as published by 0051 % the Free Software Foundation; either version 3 of the License, or 0052 % (at your option) any later version. 0053 0054 % Default values 0055 maxGap = 0.1; 0056 smooth = 2; 0057 nBins = 50; 0058 limits = [8 -2 2]; 0059 0060 % Check number of parameters 0061 if nargin < 2 | mod(length(varargin),2) ~= 0, 0062 error('Incorrect number of parameters (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0063 end 0064 0065 % Check parameter sizes 0066 if ~isdvector(spikes), 0067 error('Parameter ''spikes'' is not a vector (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0068 end 0069 if size(phases,2) ~= 2, 0070 error('Parameter ''phases'' is not a Nx2 matrix (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0071 end 0072 isradians(phases(:,2)); 0073 0074 % Parse options 0075 for i = 1:2:length(varargin), 0076 if ~ischar(varargin{i}), 0077 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).']); 0078 end 0079 switch(lower(varargin{i})), 0080 case 'smooth', 0081 smooth = varargin{i+1}; 0082 if ~isdvector(smooth,'>=0') || length(smooth) > 2, 0083 error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0084 end 0085 case 'nbins', 0086 nBins = varargin{i+1}; 0087 if ~isivector(nBins,'>0') || length(nBins) > 2, 0088 error('Incorrect value for property ''nBins'' (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0089 end 0090 case 'limits', 0091 limits = varargin{i+1}; 0092 if ~isivector(limits,'#3') || limits(1) <0 || limits(3) < limits(2), 0093 error('Incorrect value for property ''limits'' (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0094 end 0095 case 'maxgap', 0096 maxGap = varargin{i+1}; 0097 if ~isdscalar(maxGap,'>0'), 0098 error('Incorrect value for property ''maxGap'' (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).'); 0099 end 0100 otherwise, 0101 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help PhaseMango">PhaseMango</a>'' for details).']); 0102 end 0103 end 0104 0105 % Default values 0106 data.time = []; 0107 data.rate = []; 0108 data.acceleration = []; 0109 data.phase = []; 0110 map.x = []; 0111 map.y = []; 0112 map.phase = []; 0113 map.count = []; 0114 0115 % Compute spike phases 0116 if isempty(spikes), return; end 0117 spikePhases = Interpolate(phases,spikes,'trim','off','type','circular'); 0118 spikePhases = exp(j*spikePhases(:,2)); 0119 if isempty(spikePhases), return; end 0120 0121 % Determine theta cycles (start/stop) 0122 phases(:,2) = wrap(phases(:,2),1); % in [-pi,pi] 0123 start = phases(ZeroCrossings(phases),1); 0124 stop = start(2:end); 0125 start = start(1:end-1); 0126 data.time = start; 0127 0128 % Determine to which theta cycle each spike belongs 0129 [~,spikeCycles] = InIntervals(spikes,[start stop]); 0130 partialCycles = spikeCycles==0; 0131 spikeCycles(partialCycles) = []; 0132 0133 % Compute number of spikes in each cycle 0134 nSpikes = Accumulate(spikeCycles,1,size(start)); 0135 0136 % Compute average phase in each cycle 0137 spikePhases(partialCycles) = []; 0138 data.phase = angle(Accumulate(spikeCycles,spikePhases,size(start))./nSpikes); 0139 0140 phase = Accumulate(spikeCycles,spikePhases,size(start))./nSpikes; 0141 phase(nSpikes==0) = 0; 0142 kernel = ones(7,1)/7; 0143 data.phase = angle(conv(phase,kernel,'same')); 0144 0145 % Estimate spike rate and acceleration in each cycle 0146 0147 % Rate = average number of spikes over 7 surrounding cycles, which can be 0148 % efficiently computed as the convolution of nSpikes with the kernel (1/7,..,1/7) 0149 kernel = ones(7,1)/7; 0150 data.rate = conv(nSpikes,kernel,'same'); 0151 % Acceleration = slope of linear regression over 7 surrounding cycles 0152 % The general formula is b = (E[xy]-E[x]E[y]) / (E[x²]-E[x]²) 0153 % To compute this efficiently, we first note that the slope is unchanged 0154 % by translations along x the axis: we choose x in {-3..3}, so that E[x]=0. 0155 % We then compute E[x²] = 4, and end up with b = 1/4 E[xy], which can be 0156 % computed as the convolution of y with the kernel (3,2,1,0,-1,-2,-3)/7. 0157 kernel = (3:-1:-3)/7; 0158 data.acceleration = conv(nSpikes,kernel,'same')/4; 0159 0160 % Number of bins for x and y 0161 nBinsX = nBins(1); 0162 if length(nBins) == 1, 0163 nBinsY = nBinsX; 0164 nBins(2) = nBins; 0165 else 0166 nBinsY = nBins(2); 0167 end 0168 0169 % Min and max speed and acceleration 0170 xLims = [0 limits(1)]; 0171 yLims = limits(2:3); 0172 0173 % Bin rate and acceleration 0174 rate = Bin(data.rate,xLims,nBinsX); 0175 acceleration = Bin(data.acceleration,yLims,nBinsY); 0176 0177 % Compute phase map 0178 map.phase = angle(Smooth(Accumulate([rate acceleration],exp(j*data.phase),[nBinsX nBinsY]),smooth))'; 0179 map.count = Smooth(Accumulate([rate acceleration],1,[nBinsX nBinsY]),smooth)'; 0180 0181 % Bins 0182 map.x = linspace(xLims(1),xLims(2),nBinsX); 0183 map.y = linspace(yLims(1),yLims(2),nBinsY); 0184 0185 0186 % figure; 0187 % hold on; 0188 % PlotXY(phases); 0189 % PlotTicks(spikes,'size',1,'k'); 0190 % plot(start,nSpikes,'b*'); 0191 % plot(start,data.rate,'r*');