PhasePrecession - Compute spike phase precession.
Compute spike phase precession using the methods of O'Keefe and Recce (1993),
i.e. spike phase vs position, and Harris et al. (2002), i.e. spike phase vs
spike rate. Single-lap phase precession is also computed.
USAGE
[data,stats] = PhasePrecession(positions,spikes,phases,<options>)
positions linearized position <a href="matlab:help samples">samples</a> (normalized to [0..1])
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
-------------------------------------------------------------------------
'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])
'slope' search start value for slope (default = 0)
=========================================================================
NOTE
To compute only phase vs rate, positions are not required and can be left
empty, e.g. PhasePrecession([],spikes,phases).
OUTPUT
data.x position samples
data.position.ok spikes occurring at valid coordinates (logical)
data.position.t spike times (only for valid coordinates)
data.position.x x coordinate for each spike
data.position.phase spike phase for each spike (in radians)
data.position.lap lap number for each spike
data.rate.t spike times
data.rate.r spike rate for each spike
data.rate.phase spike phase for each spike (in radians)
data.rate.lap lap number for each spike
Additional statistics computed using phase vs position data:
stats.slope phase precession slope (via circular regression)
stats.intercept phase precession intercept (via circular regression)
stats.r2 coefficient of determination for circular regression
stats.p p-value for circular regression
stats.x center x for early, middle and late subfields
stats.mean mean phase for early, middle and late subfields
stats.var phase variance for early, middle and late subfields
stats.std phase standard deviation for early, middle and late subfields
stats.conf 95% confidence intervals
stats.all cell array containing all phases for each subfield
(useful for population analyses)
stats.lap.slope phase precession slope (via circular regression)
stats.lap.intercept phase precession intercept (via circular regression)
stats.lap.r2 coefficient of determination for circular regression
stats.lap.p p-value for circular regression
Additional statistics computed using phase vs rate data:
stats.rate.mean mean phase for each spike rate
stats.rate.var phase variance for each spike rate
stats.rate.conf phase 95% confidence intervals for each spike rate
stats.boundaries field boundaries
SEE
See also Phase, PlotPhasePrecession.0001 function [data,stats] = PhasePrecession(positions,spikes,phases,varargin) 0002 0003 %PhasePrecession - Compute spike phase precession. 0004 % 0005 % Compute spike phase precession using the methods of O'Keefe and Recce (1993), 0006 % i.e. spike phase vs position, and Harris et al. (2002), i.e. spike phase vs 0007 % spike rate. Single-lap phase precession is also computed. 0008 % 0009 % USAGE 0010 % 0011 % [data,stats] = PhasePrecession(positions,spikes,phases,<options>) 0012 % 0013 % positions linearized position <a href="matlab:help samples">samples</a> (normalized to [0..1]) 0014 % spikes spike timestamps 0015 % phases phase <a href="matlab:help samples">samples</a> (see <a href="matlab:help Phase">Phase</a>) 0016 % <options> optional list of property-value pairs (see table below) 0017 % 0018 % ========================================================================= 0019 % Properties Values 0020 % ------------------------------------------------------------------------- 0021 % 'maxGap' time gaps between successive position samples exceeding 0022 % this threshold (e.g. undetects) will not be interpolated 0023 % (default = 100 ms) 0024 % 'boundaries' onset and offset for single-lap phase precession can be 0025 % determined either automatically based on spike count 0026 % ('count', default) or using explicit firing field 0027 % boundaries ([Xstart Xstop], where each X is in [0..1]) 0028 % 'slope' search start value for slope (default = 0) 0029 % ========================================================================= 0030 % 0031 % NOTE 0032 % 0033 % To compute only phase vs rate, positions are not required and can be left 0034 % empty, e.g. PhasePrecession([],spikes,phases). 0035 % 0036 % OUTPUT 0037 % 0038 % data.x position samples 0039 % 0040 % data.position.ok spikes occurring at valid coordinates (logical) 0041 % data.position.t spike times (only for valid coordinates) 0042 % data.position.x x coordinate for each spike 0043 % data.position.phase spike phase for each spike (in radians) 0044 % data.position.lap lap number for each spike 0045 % 0046 % data.rate.t spike times 0047 % data.rate.r spike rate for each spike 0048 % data.rate.phase spike phase for each spike (in radians) 0049 % data.rate.lap lap number for each spike 0050 % 0051 % Additional statistics computed using phase vs position data: 0052 % 0053 % stats.slope phase precession slope (via circular regression) 0054 % stats.intercept phase precession intercept (via circular regression) 0055 % stats.r2 coefficient of determination for circular regression 0056 % stats.p p-value for circular regression 0057 % stats.x center x for early, middle and late subfields 0058 % stats.mean mean phase for early, middle and late subfields 0059 % stats.var phase variance for early, middle and late subfields 0060 % stats.std phase standard deviation for early, middle and late subfields 0061 % stats.conf 95% confidence intervals 0062 % stats.all cell array containing all phases for each subfield 0063 % (useful for population analyses) 0064 % 0065 % stats.lap.slope phase precession slope (via circular regression) 0066 % stats.lap.intercept phase precession intercept (via circular regression) 0067 % stats.lap.r2 coefficient of determination for circular regression 0068 % stats.lap.p p-value for circular regression 0069 % 0070 % Additional statistics computed using phase vs rate data: 0071 % 0072 % stats.rate.mean mean phase for each spike rate 0073 % stats.rate.var phase variance for each spike rate 0074 % stats.rate.conf phase 95% confidence intervals for each spike rate 0075 % 0076 % stats.boundaries field boundaries 0077 % 0078 % SEE 0079 % 0080 % See also Phase, PlotPhasePrecession. 0081 0082 % Copyright (C) 2004-2012 by Michaël Zugaro 0083 % 0084 % This program is free software; you can redistribute it and/or modify 0085 % it under the terms of the GNU General Public License as published by 0086 % the Free Software Foundation; either version 3 of the License, or 0087 % (at your option) any later version. 0088 0089 % Default values 0090 maxGap = 0.1; 0091 boundaries = 'count'; 0092 slope = 0; 0093 0094 % Check number of parameters 0095 if nargin < 3 | mod(length(varargin),2) ~= 0, 0096 error('Incorrect number of parameters (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).'); 0097 end 0098 0099 % Check parameter sizes 0100 if ~isempty(positions) && size(positions,2) < 2, 0101 error('Parameter ''positions'' should have at least 2 columns (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).'); 0102 end 0103 if ~isdvector(spikes), 0104 error('Parameter ''spikes'' is not a vector (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).'); 0105 end 0106 if size(phases,2) ~= 2, 0107 error('Parameter ''phases'' is not a Nx2 matrix (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).'); 0108 end 0109 isradians(phases(:,2)); 0110 0111 % Parse options 0112 for i = 1:2:length(varargin), 0113 if ~ischar(varargin{i}), 0114 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).']); 0115 end 0116 switch(lower(varargin{i})), 0117 case 'boundaries', 0118 boundaries = lower(varargin{i+1}); 0119 if ~isdvector(boundaries,'#2','>=0','<=1') && ~isastring(boundaries,'count'), 0120 error('Incorrect value for property ''boundaries'' (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).'); 0121 end 0122 case 'maxgap', 0123 maxGap = varargin{i+1}; 0124 if ~isdscalar(maxGap,'>=0'), 0125 error('Incorrect value for property ''maxGap'' (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).'); 0126 end 0127 case 'slope', 0128 slope = varargin{i+1}; 0129 if ~isdscalar(slope), 0130 error('Incorrect value for property ''slope'' (type ''help <a href="matlab:help CircularRegression">PhasePrecession</a>'' for details).'); 0131 end 0132 otherwise, 0133 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help PhasePrecession">PhasePrecession</a>'' for details).']); 0134 end 0135 end 0136 0137 % Default values 0138 data.x = positions; 0139 data.position.t = []; 0140 data.position.x = []; 0141 data.position.phase = []; 0142 data.position.lap = []; 0143 data.rate.t = []; 0144 data.rate.r = []; 0145 data.rate.phase = []; 0146 data.rate.lap = []; 0147 stats.slope = nan; 0148 stats.intercept = nan; 0149 stats.r2 = nan; 0150 stats.p = nan; 0151 stats.x = nan(1,3); 0152 stats.mean = nan(1,3); 0153 stats.var = nan(1,3); 0154 stats.std = nan(1,3); 0155 stats.conf = nan(2,3); 0156 stats.all = {nan nan nan}; 0157 stats.rate.mean = nan(1,16); 0158 stats.rate.conf = nan(2,16); 0159 stats.boundaries = []; 0160 stats.lap.slope = []; 0161 stats.lap.intercept = []; 0162 stats.lap.r2 = []; 0163 stats.lap.p = []; 0164 0165 % Compute spike phases 0166 if isempty(spikes), return; end 0167 spikePhases = Interpolate(phases,spikes,'trim','off','type','circular'); 0168 if isempty(spikePhases), return; end 0169 0170 % Interpolate positions at spike times 0171 if ~isempty(positions), 0172 % Make sure positions are normalized 0173 if max(positions(:,2)) > 1 || min(positions(:,2)) < 0, 0174 positions(:,2) = ZeroToOne(positions(:,2)); 0175 warning('Parameter ''positions'' should contain values in [0 1]. The data will now be transformed accordingly.'); 0176 end 0177 [x,ignored] = Interpolate(positions,spikes,'trim','off','maxGap',maxGap); 0178 data.position.ok = ~ignored; 0179 data.position.t = spikes(~ignored); 0180 data.position.x = x(:,2); 0181 data.position.phase = spikePhases(~ignored,2); 0182 end 0183 0184 % Compute unwrapped spike phases 0185 unwrapped = [phases(:,1) unwrap(phases(:,2))]; 0186 [spikeUnwrappedPhases,ignored] = Interpolate(unwrapped,spikes,'trim','off'); 0187 % Get beginning and end unwrapped phases of two cycles surrounding each spike 0188 startUnwrappedPhases = spikeUnwrappedPhases(:,2)-2*pi; 0189 stopUnwrappedPhases = spikeUnwrappedPhases(:,2)+2*pi; 0190 0191 %%%% 0192 % figure;PlotXY(unwrapped);hold on; 0193 % PlotHVLines(spikes,'v','k');PlotXY(spikeUnwrappedPhases,'r+');PlotHVLines(startUnwrappedPhases,'h','g');PlotHVLines(stopUnwrappedPhases,'h','r'); 0194 % figure;PlotXY(phases);hold on; 0195 % PlotTicks(spikes,pi,'k'); 0196 %%%% 0197 0198 % Compute instantaneous rate at spike times 0199 data.rate.r = CountInIntervals(spikeUnwrappedPhases(:,2),[startUnwrappedPhases stopUnwrappedPhases]); 0200 data.rate.t = spikes; 0201 data.rate.phase = spikePhases(:,2); 0202 0203 0204 % Determine lap number for each spike (for single-lap phase precession) 0205 0206 if isdvector(boundaries,'#2') && ~isempty(positions), 0207 0208 stats.boundaries = boundaries; 0209 0210 % Firing curve boundaries 0211 fieldStart = boundaries(1); 0212 fieldStop = boundaries(2); 0213 fieldSize = abs(fieldStart-fieldStop); 0214 0215 % Let x1 be current position, x0 be previous position and x2 be next position 0216 x0 = [positions(1,2);positions(1:end-1,2)]; 0217 x1 = positions(:,2); 0218 x2 = [positions(2:end,2);positions(end,2)]; 0219 0220 % Determine when the animal enters/leaves the field: 0221 % 1) the animal enters the field when x1 is inside and x0 is not inside and near x1 0222 % 2) the animal leaves the field when x1 is inside and x2 is not inside and near x1 0223 % There are two additional cases: 0224 % 1) the animal enters the field in the first lap when x1 is inside and close to the entrance 0225 % 1) the animal leaves the field in the last lap when x1 is inside and close to the exit 0226 if fieldStart < fieldStop, 0227 x0_inside = x0 >= fieldStart & x0 <= fieldStop; 0228 x1_inside = x1 >= fieldStart & x1 <= fieldStop; 0229 x2_inside = x2 >= fieldStart & x2 <= fieldStop; 0230 x0_near_x1 = abs(x0-x1)<0.1*fieldSize; 0231 x2_near_x1 = abs(x2-x1)<0.1*fieldSize; 0232 else 0233 x0_inside = x0 >= fieldStart | x0 <= fieldStop; 0234 x1_inside = x1 >= fieldStart | x1 <= fieldStop; 0235 x2_inside = x2 >= fieldStart | x2 <= fieldStop; 0236 x0_near_x1 = abs(x0-x1)<0.1*fieldSize | 1-abs(x0-x1)<0.1*fieldSize; 0237 x2_near_x1 = abs(x2-x1)<0.1*fieldSize | 1-abs(x2-x1)<0.1*fieldSize; 0238 end 0239 start = x1_inside & ~(x0_inside & x0_near_x1); 0240 start(1) = x1_inside(1) & abs(x1(1)-fieldStart)<0.1*fieldSize; 0241 stop = x1_inside & ~(x2_inside & x2_near_x1); 0242 stop(end) = x1_inside(end) & abs(x1(end)-fieldStop)<0.1*fieldSize; 0243 0244 if any(start) && any(stop), 0245 startT = positions(start,1); 0246 stopT = positions(stop,1); 0247 % Get rid of duplicate events (e.g. the rat enters twice before it exits) 0248 s = [startT ones(size(startT));stopT 2*ones(size(stopT))]; 0249 s = sortrows(s); 0250 ds = [1;diff(s(:,2))]; 0251 s = s(ds~=0,:); 0252 startT = s(s(:,2)==1); 0253 stopT = s(s(:,2)==2); 0254 % Drop incomplete traversals of the field 0255 if startT(1) > stopT(1), stopT(1) = []; end 0256 if ~isempty(stopT), 0257 if startT(end) > stopT(end), startT(end) = []; end 0258 traversals = [startT stopT]; 0259 0260 % Show laps (debugging) 0261 %figure;hold on; 0262 %nn = positions(:,1); 0263 %plot(nn,x1); 0264 %PlotHVLines([fieldStart fieldStop],'h',':k'); 0265 %PlotHVLines(startT,'v','g'); 0266 %PlotHVLines(stopT,'v','r'); 0267 0268 % Determine lap number 0269 [~,lap] = InIntervals(data.position.t,traversals); 0270 data.position.lap = lap; 0271 [~,lap] = InIntervals(data.rate.t,traversals); 0272 data.rate.lap = lap; 0273 end 0274 end 0275 0276 elseif strcmp(boundaries,'count'), 0277 0278 % Get beginning and end unwrapped phases of eight cycles before each spike 0279 startUnwrappedPhases = spikeUnwrappedPhases(:,2)-8*2*pi; 0280 stopUnwrappedPhases = spikeUnwrappedPhases(:,2); 0281 rateBefore = CountInIntervals(spikeUnwrappedPhases,[startUnwrappedPhases stopUnwrappedPhases]); 0282 % Get beginning and end unwrapped phases of eight cycles after each spike 0283 startUnwrappedPhases = spikeUnwrappedPhases(:,2); 0284 stopUnwrappedPhases = spikeUnwrappedPhases(:,2)+8*2*pi; 0285 rateAfter = CountInIntervals(spikeUnwrappedPhases,[startUnwrappedPhases stopUnwrappedPhases]); 0286 % Find periods of minimal (resp. maximal) firing rates 0287 minimalBefore = rateBefore <= 4; 0288 maximalAfter = rateAfter >= 16; 0289 maximalBefore = rateBefore >= 16; 0290 minimalAfter = rateAfter <= 4; 0291 start = minimalBefore & maximalAfter; 0292 stop = maximalBefore & minimalAfter; 0293 % Drop incomplete traversals of the field 0294 if any(start) && any(stop), 0295 if start(1) >= stop(1), stop(1) = []; end 0296 if start(end) >= stop(end), start(end) = []; end 0297 traversals = [start stop]; 0298 % Determine lap number 0299 [~,lap] = InIntervals(data.position.t,traversals); 0300 data.position.lap = lap; 0301 [~,lap] = InIntervals(data.rate.t,traversals); 0302 data.rate.lap = lap; 0303 end 0304 0305 end 0306 0307 % Statistics: regression + circular means and variances for early, middle and late subfields 0308 if ~isempty(positions) && nargout >= 2, 0309 % For circular/circular regression, shift the field leftward along x (if necessary) 0310 % to treat as circular-linear regression - as a result, the intercept corresponds to 0311 % the right portion of the field 0312 x = data.position.x; 0313 dx = 0; 0314 if isdvector(boundaries,'#2'), 0315 if fieldStart > fieldStop, 0316 dx = fieldStop; % shift by this amount 0317 left = x < fieldStart; 0318 x(left) = x(left) + 1 - dx; 0319 x(~left) = x(~left) - dx; 0320 end 0321 end 0322 % Circular regression for average data 0323 ok = ~isnan(x); 0324 [beta,r2,p] = CircularRegression(x(ok),data.position.phase(ok),'slope',slope); 0325 stats.slope = beta(1); 0326 stats.intercept = beta(2) + stats.slope*dx; % Correction for circular/circular regression (see above) 0327 stats.r2 = r2; 0328 stats.p = p; 0329 % Circular regression for each lap separately 0330 for lap = 1:max(data.rate.lap), 0331 thisLap = data.position.lap == lap; 0332 [beta,r2,p] = CircularRegression(x(thisLap&ok),data.position.phase(thisLap&ok),'slope',stats.slope); 0333 stats.lap.slope(lap,1) = beta(1); 0334 stats.lap.intercept(lap,1) = beta(2) + stats.lap.slope(lap,1)*dx; 0335 stats.lap.r2(lap,1) = r2; 0336 stats.lap.p(lap,1) = p; 0337 end 0338 % Mean phase for each subfield 0339 x0 = min(data.position.x); 0340 x1 = max(data.position.x); 0341 dx = x1-x0; 0342 for i = 1:3, 0343 ok = data.position.x > x0+dx/3*(i-1) & data.position.x < x0+dx/3*i; 0344 if sum(ok) == 0, continue; end 0345 stats.all{i} = data.position.phase(ok); 0346 stats.x(i) = x0+1.5*dx/3*(i-1); 0347 [stats.var(i),stats.std(i)] = CircularVariance(data.position.phase(ok)); 0348 [stats.mean(i),stats.conf(:,i)] = CircularConfidenceIntervals(data.position.phase(ok)); 0349 end 0350 end 0351 0352 % Statistics: phase vs rate 0353 if nargout >= 2, 0354 for rate = 1:max(data.rate.r), 0355 ok = data.rate.r==rate; 0356 if any(ok), 0357 [stats.rate.mean(rate),stats.rate.conf(:,rate)] = CircularConfidenceIntervals(data.rate.phase(ok)); 0358 else 0359 stats.rate.mean(rate) = NaN; 0360 stats.rate.conf(1:2,rate) = NaN; 0361 end 0362 end 0363 end