TestRemapping

TestRemapping - Test if firing fields remap (or shift) between two conditions.

  To test for random remapping between a control and a test condition, we
  first estimate for each place cell the absolute field shift between the two
  conditions, as the mode of the spatial cross-correlogram between the
  respective firing fields. We then divide this value by the average field
  size across conditions, yielding a relative shift. This measures by how much
  the field moves between the two conditions, in proportion to the field size.
  To account for baseline variability, we also estimate the relative field
  shift between two repetitions of the control condition, and subtract it from
  the relative shift in the test condition. This (unsigned) corrected shift is
  averaged over all cells.

  To test if this is significantly different from random remapping, we generate
  a distribution of corrected relative shifts under the null hypothesis of
  random remapping. This is done by randomly permuting the test firing fields
  between the cells (bootstrap).

  The null hypothesis is rejected if the probability of the observed shift is
  lower than the critical value.

  In addition, we compute the (1-alpha) confidence interval of the (signed)
  corrected shifts. These can be used to test for systematic (expected) field
  shifts.

  Current implementation only tests 1D environments.

  Note that a repeated control condition is not required: one may provide the
  same control condition twice, in which case the correction will be zero.

  USAGE

    [h,p,cr,ca] = TestRemapping(control,repeat,test,<options>)

    control        firing fields (rates) in control condition (MxN: M fields, N bins)
    repeat         firing fields in repeated control condition
    test           firing fields in test condition
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'type'        'linear' for linear tracks (default), 'circular' otherwise
     'alpha'       significance level (default = 0.05)
     'iterations'  number of iterations for bootstrap (default = 150)
     'show'        plot results (default = 'off')
    =========================================================================

  OUTPUT

    h              0 if H0 (random remapping) cannot be rejected, 1 otherwise
    p              p-value of bootstrap test
    cr             confidence interval for relative shifts
    ca             confidence interval for absolute shifts

0001 function [h,p,cr,ca] = TestRemapping(control,repeat,test,varargin)
0002 
0003 %TestRemapping - Test if firing fields remap (or shift) between two conditions.
0004 %
0005 %  To test for random remapping between a control and a test condition, we
0006 %  first estimate for each place cell the absolute field shift between the two
0007 %  conditions, as the mode of the spatial cross-correlogram between the
0008 %  respective firing fields. We then divide this value by the average field
0009 %  size across conditions, yielding a relative shift. This measures by how much
0010 %  the field moves between the two conditions, in proportion to the field size.
0011 %  To account for baseline variability, we also estimate the relative field
0012 %  shift between two repetitions of the control condition, and subtract it from
0013 %  the relative shift in the test condition. This (unsigned) corrected shift is
0014 %  averaged over all cells.
0015 %
0016 %  To test if this is significantly different from random remapping, we generate
0017 %  a distribution of corrected relative shifts under the null hypothesis of
0018 %  random remapping. This is done by randomly permuting the test firing fields
0019 %  between the cells (bootstrap).
0020 %
0021 %  The null hypothesis is rejected if the probability of the observed shift is
0022 %  lower than the critical value.
0023 %
0024 %  In addition, we compute the (1-alpha) confidence interval of the (signed)
0025 %  corrected shifts. These can be used to test for systematic (expected) field
0026 %  shifts.
0027 %
0028 %  Current implementation only tests 1D environments.
0029 %
0030 %  Note that a repeated control condition is not required: one may provide the
0031 %  same control condition twice, in which case the correction will be zero.
0032 %
0033 %  USAGE
0034 %
0035 %    [h,p,cr,ca] = TestRemapping(control,repeat,test,<options>)
0036 %
0037 %    control        firing fields (rates) in control condition (MxN: M fields, N bins)
0038 %    repeat         firing fields in repeated control condition
0039 %    test           firing fields in test condition
0040 %    <options>      optional list of property-value pairs (see table below)
0041 %
0042 %    =========================================================================
0043 %     Properties    Values
0044 %    -------------------------------------------------------------------------
0045 %     'type'        'linear' for linear tracks (default), 'circular' otherwise
0046 %     'alpha'       significance level (default = 0.05)
0047 %     'iterations'  number of iterations for bootstrap (default = 150)
0048 %     'show'        plot results (default = 'off')
0049 %    =========================================================================
0050 %
0051 %  OUTPUT
0052 %
0053 %    h              0 if H0 (random remapping) cannot be rejected, 1 otherwise
0054 %    p              p-value of bootstrap test
0055 %    cr             confidence interval for relative shifts
0056 %    ca             confidence interval for absolute shifts
0057 
0058 % Copyright (C) 2012 by Michaël Zugaro
0059 %
0060 % This program is free software; you can redistribute it and/or modify
0061 % it under the terms of the GNU General Public License as published by
0062 % the Free Software Foundation; either version 3 of the License, or
0063 % (at your option) any later version.
0064 
0065 % Defaults
0066 type = 'linear';
0067 show = 'off';
0068 nBootstrap = 150;
0069 alpha = 0.05;
0070 
0071 % Check number of parameters
0072 if nargin < 3 | mod(length(varargin),2) ~= 0,
0073   error('Incorrect number of parameters (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0074 end
0075 
0076 % Check parameter sizes
0077 if ~isdmatrix(control) || ~isdmatrix(repeat) || ~isdmatrix(test),
0078     error('All firing fields should be MxN matrices (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0079 end
0080 if ~(all(size(control)==size(repeat)) && all(size(control)==size(test))),
0081     error('All firing fields should have the same sizes (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0082 end
0083 
0084 % Parse parameter list
0085 for i = 1:2:length(varargin),
0086     if ~ischar(varargin{i}),
0087         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).']);
0088     end
0089     switch(lower(varargin{i})),
0090         case 'type',
0091             type = varargin{i+1};
0092             if ~isastring(type,'linear','circular'),
0093                 error('Incorrect value for property ''type'' (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0094             end
0095         case 'alpha',
0096             alpha = varargin{i+1};
0097             if ~isdscalar(alpha,'>=0','<=1'),
0098                 error('Incorrect value for property ''alpha'' (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0099             end
0100         case 'iterations',
0101             nBootstrap = varargin{i+1};
0102             if ~isiscalar(nBootstrap,'>0'),
0103                 error('Incorrect value for property ''iterations'' (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0104             end
0105         case 'show',
0106             show = varargin{i+1};
0107             if ~isastring(show,'on','off'),
0108                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).');
0109             end
0110         otherwise,
0111             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help TestRemapping">TestRemapping</a>'' for details).']);
0112     end
0113 end
0114 [m,n] = size(control);
0115 
0116 % Compute relative field shift (= relative to size)
0117 [relativeShiftControl,absoluteShiftControl,xcControl] = FieldShift(control,repeat,'type',type);
0118 [relativeShiftTest,absoluteShiftTest,xcTest] = FieldShift(control,test,'type',type);
0119 
0120 % Mean relative shift, taking into account baseline variability
0121 meanUnsignedRelativeShift = mean(abs(relativeShiftTest));
0122 
0123 % Cumulative density function (CDF) of null distribution, estimated by bootstrap
0124 global RemappingTest_control RemappingTest_shiftControl RemappingTest_type;
0125 RemappingTest_control = control;
0126 RemappingTest_shiftControl = relativeShiftControl;
0127 RemappingTest_type = type;
0128 bRelativeShift = bootstrp(nBootstrap,@RemappingShift,control); % see RemappingShift below
0129 [bPDF,x] = hist(bRelativeShift,100);
0130 bPDF = bPDF / nBootstrap;
0131 bCDF = cumsum(bPDF);
0132 
0133 % Test for random remapping
0134 if meanUnsignedRelativeShift < min(x),
0135     p = 0;
0136 elseif meanUnsignedRelativeShift > max(x),
0137     p = 1;
0138 else
0139     p = interp1(x,bCDF,meanUnsignedRelativeShift);
0140 end
0141 h = double(p<alpha);
0142 
0143 % Confidence intervals for relative and absolute shifts (via bootstrap)
0144 s = bootstrp(nBootstrap,@mean,absoluteShiftTest);
0145 ca = [prctile(s,100*(alpha/2)) prctile(s,100*(1-(alpha/2)))];
0146 s = bootstrp(nBootstrap,@mean,relativeShiftTest);
0147 cr = [prctile(s,100*(alpha/2)) prctile(s,100*(1-(alpha/2)))];
0148 
0149 % Find peak positions (optimized code, uneasy to read)
0150 [y,x] = find(control==repmat(max(control,[],2),1,n));
0151 peakControl = Accumulate(y,x)./Accumulate(y,1)/n;
0152 [y,x] = find(test==repmat(max(test,[],2),1,n));
0153 peakTest = Accumulate(y,x)./Accumulate(y,1)/n;
0154 
0155 if strcmp(show,'on'),
0156 
0157     figure;
0158     % Plot control cross-correlograms
0159     subplot(1,2,1);hold on;
0160     PlotColorCurves(xcControl,[-1 1],0,'w--',absoluteShiftControl,'.k','gamma',2);
0161     title(['Control (N=' int2str(size(xcControl,1)) ')']);
0162     % Plot test cross-correlograms
0163     subplot(1,2,2);hold on;
0164     PlotColorCurves(xcTest,[-1 1],0,'w--',absoluteShiftTest,'.k','gamma',2);
0165     title(['Test (N=' int2str(size(xcTest,1)) ')']);
0166 
0167    figure;
0168     % Plot control cross-correlograms in ascending order
0169     [absoluteShiftControl,order] = sort(absoluteShiftControl);
0170     xcControl = xcControl(order,:);
0171     subplot(1,2,1);hold on;
0172     PlotColorCurves(xcControl,[-1 1],0,'w--',absoluteShiftControl,'.k','gamma',2);
0173     % Plot histogram
0174     SideAxes('top',0.1,'gap',0.01);
0175     [hi,xh] = hist(absoluteShiftControl,50);
0176     bar(xh,hi);
0177     PlotHVLines(mean(absoluteShiftControl),'v','r','linewidth',2);
0178     xlim([-1 1]);
0179     set(gca,'xtick',[]);
0180     title(['Control (N=' int2str(size(xcControl,1)) ')']);
0181     % Plot test cross-correlograms in ascending order
0182     [absoluteShiftTest,order] = sort(absoluteShiftTest);
0183     xcTest = xcTest(order,:);
0184     subplot(1,2,2);hold on;
0185     PlotColorCurves(xcTest,[-1 1],0,'w--',absoluteShiftTest,'.k','gamma',2);
0186     % Plot histogram
0187     SideAxes('top',0.1,'gap',0.01);
0188     [hi,xh] = hist(absoluteShiftTest,50);
0189     bar(xh,hi);
0190     PlotHVLines(mean(absoluteShiftTest),'v','r','linewidth',2);
0191     xlim([-1 1]);
0192     set(gca,'xtick',[]);
0193     title(['Test (N=' int2str(size(xcTest,1)) ')']);
0194 
0195     figure;hold on;
0196     % Plot relative shift histogram
0197     [hi,xh] = hist(relativeShiftTest,50);
0198     bar(xh,hi/nBootstrap);
0199     PlotHVLines(cr,'v','k--');
0200     PlotHVLines(mean(s),'v','r','linewidth',2);
0201     xlabel('Relative Shift');
0202     ylabel('Probability');
0203     title(['N = ' num2str(nBootstrap)]);
0204     xLim = xlim;
0205     yLim = ylim;
0206     [~,pval] = ttest(relativeShiftTest);
0207     prop1 = sum(abs(relativeShiftTest)<=1)/length(relativeShiftTest);
0208     prop2 = sum(abs(relativeShiftTest)<=0.5)/length(relativeShiftTest);
0209     m = mean(relativeShiftTest);
0210     s = sem(relativeShiftTest);
0211     l = length(relativeShiftTest);
0212     text(xLim(1)+0.1*abs(xLim(1)),yLim(2)*0.9,{mean2str(m,s,l),[num2str((1-alpha)*100) '% Confidence Interval [' sprintf('%.3f',cr(1)) ',' num2str(cr(2)) ']'],['t-test p = ' sprintf('%.3f',pval)],[sprintf('%.3f',prop1*100) '% |x|<=1'],[sprintf('%.3f',prop2*100) '% |x|<=0.5']});
0213 
0214     figure;
0215     % Plot test vs control shifts (density plot)
0216     subplot(2,2,1);hold on;
0217     nBins = 500;
0218     range = max([relativeShiftControl;relativeShiftTest])*[-1 1];
0219     sc = Bin(relativeShiftControl,range,nBins);
0220     st = Bin(relativeShiftTest,range,nBins);
0221     z = Smooth(Accumulate([sc st],1,[nBins nBins]),10)'/n;
0222     PlotColorMap(z,'x',linspace(range(1),range(2),n),'y',linspace(range(1),range(2),n));
0223     xlabel('Control Shift (proportion of field size)');
0224     ylabel('Test Shift (proportion of field size)');
0225     clim([0 0.3*max(z(:))]);
0226     % Plot field modes in control vs test (density plot)
0227     subplot(2,2,2);hold on;
0228     PlotRepeat(peakControl,peakTest,'xxyy','k.');
0229     xlabel('Control Peak Position');
0230     ylabel('Test Peak Position');
0231     [beta,R2,pSlope] = CircularRegression(peakControl*2*pi,peakTest*2*pi,'slope',1);
0232     title(['slope = ' num2str(beta(1)) ' (p = ' num2str(pSlope) ')']);
0233     PlotSlope(1,1,beta(1),2,'r','LineWidth',2);
0234     % Plot shift distribution under H0 (bootstrap)
0235     subplot(2,2,[3 4]);hold on;
0236     [hi,x] = hist(bRelativeShift,30);
0237     b = bar(x,hi/nBootstrap);
0238     hold on;
0239     PlotHVLines(meanUnsignedRelativeShift,'v','r');
0240     Y = ylim; Y = Y(2);
0241     text(meanUnsignedRelativeShift*1.25,Y*0.85,['p = ' sprintf('%.3f',p)]);
0242     xlim([0 1]);
0243     title(['N = ' num2str(nBootstrap)]);
0244     xlabel('Mean Unsigned Shift (proportion of field size)');
0245     ylabel('Probability');
0246     
0247     figure;
0248     % Plot shift as a function of peak position
0249     subplot(2,2,1);hold on;
0250     pc = Bin(peakControl,[0 1],10);
0251     [m,stdev] = Accumulate(pc,relativeShiftTest);
0252     s = stdev/sqrt(length(m));
0253     %errorbar(0.05+(0:0.1:0.9),m,ci(:,1)-m,ci(:,2)-m,'k','LineStyle','none');
0254     errorbar(0.05+(0:0.1:0.9),m,s,s,'k','LineStyle','none');
0255     bar(0.05+(0:0.1:0.9),m);
0256     plot(peakControl,relativeShiftTest,'.r');
0257     xlabel('Peak Position');
0258     ylabel('Shift (proportion of field size)');
0259     pa = anova1(relativeShiftTest,pc,'off');
0260     yLim = ylim;
0261     text(mean(xlim),yLim(2)*0.9,['p = ' num2str(pa)]);
0262     subplot(2,2,3);hold on;
0263     [m,stdev] = Accumulate(pc,abs(relativeShiftTest));
0264     s = stdev/sqrt(length(m));
0265     errorbar(0.05+(0:0.1:0.9),m,s,s,'k','LineStyle','none');
0266     bar(0.05+(0:0.1:0.9),m);
0267     plot(peakControl,relativeShiftTest,'.r');
0268     xlabel('Peak Position');
0269     ylabel('Unsigned Shift (proportion of field size)');
0270     pa = anova1(abs(relativeShiftTest),pc,'off');
0271     yLim = ylim;
0272     text(mean(xlim),yLim(2)*0.9,['p = ' num2str(pa)]);
0273     subplot(2,2,2);hold on;
0274     [m,stdev] = Accumulate(pc,absoluteShiftTest);
0275     s = stdev/sqrt(length(m));
0276     errorbar(0.05+(0:0.1:0.9),m,s,s,'k','LineStyle','none');
0277     bar(0.05+(0:0.1:0.9),m);
0278     plot(peakControl,absoluteShiftTest,'.r');
0279     xlabel('Peak Position');
0280     ylabel('Shift (proportion of field size)');
0281     pa = anova1(absoluteShiftTest,pc,'off');
0282     yLim = ylim;
0283     text(mean(xlim),yLim(2)*0.9,['p = ' num2str(pa)]);
0284     subplot(2,2,4);hold on;
0285     [m,stdev] = Accumulate(pc,abs(absoluteShiftTest));
0286     s = stdev/sqrt(length(m));
0287     errorbar(0.05+(0:0.1:0.9),m,s,s,'k','LineStyle','none');
0288     bar(0.05+(0:0.1:0.9),m);
0289     plot(peakControl,absoluteShiftTest,'.r');
0290     xlabel('Peak Position');
0291     ylabel('Unsigned Shift (proportion of field size)');
0292     pa = anova1(abs(absoluteShiftTest),pc,'off');
0293     yLim = ylim;
0294     text(mean(xlim),yLim(2)*0.9,['p = ' num2str(pa)]);
0295 
0296 end
0297 
0298 % ------------------------------- Helper function -------------------------------
0299 
0300 function meanUnsignedRelativeShift = RemappingShift(test)
0301 
0302 global RemappingTest_control RemappingTest_shiftControl RemappingTest_type;
0303 control = RemappingTest_control;
0304 relativeShiftControl = RemappingTest_shiftControl;
0305 type = RemappingTest_type;
0306 
0307 % Compute field shift as the distance relative to size
0308 relativeShiftTest = FieldShift(control,test,'type',type);
0309 
0310 % Shift, taking into account 'intrinsic' variability
0311 meanUnsignedRelativeShift = mean(abs(relativeShiftTest));