Description of FieldPSP

0001 function stats = FieldPSP(stims,lfp,varargin)
0002 
0003 %FieldPSP - Measure field EPSPs (waveforms, slope, amplitude) across time.
0004 %
0005 %  Slope is determined as the maximum slope during the time window around stimulation.
0006 %  Amplitude is determined as the maximum amplitude during the same time window.
0007 %
0008 %  USAGE
0009 %
0010 %    stats = FieldPSP(stims,lfp,<options>)
0011 %
0012 %    stims          stimulation timestamps
0013 %    lfp            local field potential <a href="matlab:help samples">samples</a>
0014 %    <options>      optional list of property-value pairs (see table below)
0015 %
0016 %    =========================================================================
0017 %     Properties    Values
0018 %    -------------------------------------------------------------------------
0019 %     'direction'   direction of the field EPSP, either 'up' (recording
0020 %                   from str. pyr., default) or 'down' (recording from str.
0021 %                   radiatum)
0022 %     'durations'   durations before and after synchronizing events for each
0023 %                   trial (in s) (default = [-0.01 0.02])
0024 %     'slope'       where max slope should be sought (default = [0.001 0.005])
0025 %     'amplitude'   where max amplitude should be sought
0026 %                   (default = [0.005 0.010])
0027 %     'mode'        compute slope and amplitude running averages using a
0028 %                   fixed time window ('time'), or a window with a fixed
0029 %                   number of stimulations ('count', default)
0030 %     'window'      window length (in s or in counts) for running averages
0031 %                   (default = 100)
0032 %     'show'        either 'on' or 'off' (default)
0033 %    =========================================================================
0034 %
0035 %  OUTPUT
0036 %
0037 %    stats.amplitude.t         times of max amplitude (relative to stimulation)
0038 %    stats.amplitude.v         values of max amplitudes
0039 %    stats.slope.t             times of max slopes (relative to stimulation)
0040 %    stats.slope.v             values of max slopes
0041 %    stats.psp.t               times of fPSP samples (relative to stimulation)
0042 %    stats.psp.v               fPSP waveforms, one per stimulation
0043 %
0044 
0045 % Copyright (C) 2010-2012 by Michaël Zugaro
0046 %
0047 % This program is free software; you can redistribute it and/or modify
0048 % it under the terms of the GNU General Public License as published by
0049 % the Free Software Foundation; either version 3 of the License, or
0050 % (at your option) any later version.
0051 
0052 % Default values
0053 parent = [];
0054 durations = [-0.01 0.02];
0055 mode = 'count';
0056 window = 100;
0057 direction = 'up';
0058 stats = [];
0059 show = 'off';
0060 slopeWindow = [0.001 0.005];
0061 amplitudeWindow = [0.005 0.010];
0062 
0063 % Check number of parameters
0064 if nargin < 2 | mod(length(varargin),2) ~= 0,
0065   error('Incorrect number of parameters (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0066 end
0067 
0068 % Check parameter sizes
0069 if ~isdvector(stims),
0070     error('Parameter ''stims'' is not a vector (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0071 end
0072 if ~isdmatrix(lfp),
0073     error('Parameter ''lfp'' is not a matrix (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0074 end
0075 
0076 % Parse parameter list
0077 for i = 1:2:length(varargin),
0078     if ~ischar(varargin{i}),
0079         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).']);
0080     end
0081     switch(lower(varargin{i})),
0082         case 'direction',
0083             direction = lower(varargin{i+1});
0084             if ~isastring(direction,'up','down'),
0085                 error('Incorrect value for property ''direction'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0086             end
0087         case 'durations',
0088             durations = varargin{i+1};
0089             if ~isdvector(durations,'#2','<'),
0090                 error('Incorrect value for property ''durations'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0091             end
0092         case 'slope',
0093             slopeWindow = varargin{i+1};
0094             if ~isdvector(slopeWindow,'#2','<'),
0095                 error('Incorrect value for property ''slope'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0096             end
0097         case 'amplitude',
0098             amplitudeWindow = varargin{i+1};
0099             if ~isdvector(amplitudeWindow,'#2','<'),
0100                 error('Incorrect value for property ''amplitude'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0101             end
0102         case 'mode',
0103             mode = lower(varargin{i+1});
0104             if ~isastring(mode,'time','count'),
0105                 error('Incorrect value for property ''mode'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0106             end
0107         case 'parent',
0108             parent = varargin{i+1};
0109             if ~ishandle(parent),
0110                 error('Incorrect value for property ''parent'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0111             end
0112         case 'window',
0113             window = varargin{i+1};
0114             if ~isdscalar(window,'>0'),
0115                 error('Incorrect value for property ''window'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0116             end
0117         case 'show',
0118             show = lower(varargin{i+1});
0119             if ~isastring(show,'on','off'),
0120                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).');
0121             end
0122         otherwise,
0123             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help FieldPSP">FieldPSP</a>'' for details).']);
0124     end
0125 end
0126 
0127 up = strcmp(direction,'up');
0128 show = strcmp(show,'on');
0129 
0130 % Some analyses depend on the method (fixed time vs fixed count), but others do not
0131 % To simplify the code, we will use two variables:
0132 %  - stims0 uses time, whatever the method
0133 %  - stims  uses either time or count, depending on the method
0134 stims0 = stims;
0135 if strcmp(mode,'count'),
0136     stims = (1:length(stims))';
0137 end
0138 
0139 % Constants
0140 nStims = length(stims);
0141 height = max(lfp(:,2))-min(lfp(:,2));
0142 nSlices = 12;
0143 start = stims(1);
0144 stop = stims(end);
0145 window = (stop-start)/nSlices;
0146 slices = (0:nSlices-1)*window+start;
0147 nStims = length(stims);
0148 frequency = 1/median(diff(lfp(:,1)));
0149 smooth = frequency/1250;
0150 nStimsPlottedPerSlice = 5;
0151 nStimsPerSlice = nStims/nSlices;
0152 nSkip = floor(nStims/(nSlices*(nStimsPlottedPerSlice-1))); % plot one stim every n
0153 colors = Bright(nStimsPlottedPerSlice);
0154 
0155 % Compute slope and amplitude for each stimulation
0156 
0157 % Resynchronize all fPSPs (use only 15ms after stims)
0158 [sync,index] = Sync(lfp,stims0,'durations',durations);
0159 
0160 % Preallocate large matrix to store all fPSPs
0161 n = sum(index==1);
0162 stats.psp.v = nan(nStims,n);
0163 stats.psp.t = linspace(durations(1),durations(2),n);
0164 % Compute slope and amplitude for each stimulation
0165 if show, figure; end
0166 for i = 1:nStims,
0167     % Select the appropriate data, store in fPSP matrix, smooth and differentiate
0168     s = sync(index==i,:);
0169     stats.psp.v(i,:) = interp1(s(:,1),s(:,2),stats.psp.t);
0170     s(:,2) = Smooth(s(:,2),smooth);
0171     % Compute amplitude as the first local maximum (or minimum) in the appropriate time interval
0172     in = InIntervals(s(:,1),amplitudeWindow);
0173     if up,
0174         extrema = find(IsExtremum(s)&in);
0175     else
0176         extrema = find(IsExtremum(s,'mode','minima')&in);
0177     end
0178     if isempty(extrema),
0179         amplitude(i) = nan;
0180         tAmplitude(i) = nan;
0181         slope(i) = nan;
0182         tSlope(i) = nan;
0183         continue;
0184     end
0185     extremum = extrema(1);
0186     tAmplitude(i) = s(extremum,1);
0187     amplitude(i) = s(extremum,2);
0188     % Compute slope as the maximum slope in the appropriate time interval
0189     % (set all values outside this time interval to nan to make sure the max will be in the interval)
0190     ds = Diff(s);
0191     ds(~InIntervals(ds(:,1),slopeWindow),2) = nan;
0192     if up,
0193         [slope(i),j] = nanmax(ds(:,2));
0194     else
0195         [slope(i),j] = nanmin(ds(:,2));
0196     end
0197     tSlope(i) = ds(j,1);
0198     % Plot slope and amplitude
0199     slice = floor(i/nStimsPerSlice);
0200     stimInSlice = floor(i-1-slice*nStimsPerSlice);
0201     if show && rem(stimInSlice,nSkip) == 0 && stimInSlice >= 0,
0202         k = ceil(nSlices*i/nStims);
0203         l = floor(stimInSlice/nSkip)+1;
0204         if l <= nStimsPlottedPerSlice,
0205             SquareSubplot(nSlices,k);
0206             xlabel(num2str(slices(k)));
0207             hold on;
0208             PlotXY(s,'color',colors(l,:),'linewidth',3);
0209             PlotSlope(tAmplitude(i),amplitude(i),0,0.005,'r');
0210             PlotSlope(tSlope(i),s(j,2),slope(i),0.005,'k');
0211         end
0212     end
0213 end
0214 
0215 % Make sure all subplots use the same scale
0216 if show, AdjustAxes('y','uniform'); end
0217 
0218 stats.amplitude.t = tAmplitude;
0219 stats.amplitude.v = amplitude;
0220 stats.slope.t = tSlope;
0221 stats.slope.v = slope;
0222 
0223 % Additional figures: mean fPSPs, slopes and amplitudes (scatterplots + running averages), CVs, ISIs and short time CCGs
0224 
0225 if show,
0226     % Mean fPSPs
0227     figure;
0228     for i = 1:nSlices,
0229         if strcmp(mode,'count'),
0230             start = floor((i-1)*nStimsPerSlice)+1;
0231             stop = floor(i*nStimsPerSlice);
0232             ok = start:stop;
0233         else
0234             ok = find(InIntervals(stims,[0 window]+start+(i-1)*window));
0235         end
0236         m = mean(stats.psp.v(ok,:));
0237         e = sem(stats.psp.v(ok,:));
0238         SquareSubplot(nSlices,i);
0239         xlabel(num2str(slices(i)));
0240         hold on;
0241         PlotMean(stats.psp.t,m,m-e/2,m+e/2,':');
0242     end
0243     % Make sure all subplots use the same scale
0244     AdjustAxes('y','uniform');
0245 
0246     figure;
0247     % 'Rectify' slope and amplitude
0248     if ~up,
0249         amplitude = -amplitude;
0250         slope = -slope;
0251     end
0252     if strcmp(mode,'time'),
0253         xLabel = 'Time (s)';
0254     else
0255         xLabel = 'Count';
0256     end
0257     % Slope (scatterplot)
0258     subplot(2,3,1);
0259     hold on;
0260     plot(stims,slope,'.','markersize',1);
0261     % Slope (running average)
0262     [ts,ms,es] = RunningAverage(stims,slope,'window',window,'overlap',0.5*window);
0263     PlotMean(ts,ms,es(:,1),es(:,2),':','k');
0264     PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
0265     xlim([stims(1,1) stims(end,1)]);
0266     ylabel('Slope (a.u.)');
0267     set(gca,'xtick',[]);
0268     % Slope time (scatterplot)
0269     subplot(2,3,2);
0270     hold on;
0271     plot(stims,tSlope*1000,'.','markersize',1);
0272     % Slope time (running average)
0273     [ts,ms,es] = RunningAverage(stims,tSlope*1000,'window',window,'overlap',0.5*window);
0274     PlotMean(ts,ms,es(:,1),es(:,2),':','k');
0275     PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
0276     xlim([stims(1,1) stims(end,1)]);
0277     ylabel('Slope time (ms)');
0278     set(gca,'xtick',[]);
0279 
0280     % Amplitude (scatter)
0281     subplot(2,3,4);
0282     hold on;
0283     plot(stims,amplitude,'.','markersize',1);
0284     % Amplitude (running average)
0285     [ta,ma,ea] = RunningAverage(stims,amplitude,'window',window,'overlap',0.5*window);
0286     PlotMean(ta,ma,ea(:,1),ea(:,2),':','k');
0287     PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
0288     xlim([stims(1,1) stims(end,1)]);
0289     xlabel(xLabel);
0290     ylabel('Amplitude (a.u.)');
0291     % Amplitude time (scatter)
0292     subplot(2,3,5);
0293     hold on;
0294     plot(stims,tAmplitude*1000,'.','markersize',1);
0295     % Amplitude time (running average)
0296     [ta,ma,ea] = RunningAverage(stims,tAmplitude*1000,'window',window,'overlap',0.5*window);
0297     xlim([stims(1,1) stims(end,1)]);
0298     PlotMean(ta,ma,ea(:,1),ea(:,2),':','k');
0299     PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
0300     ylabel('Amplitude time (ms)');
0301     xlabel(xLabel);
0302 
0303     % ISIs
0304     subplot(2,3,3);
0305     ds = diff(stims0(:,1));
0306     x = 0:0.1:10;
0307     h = hist(ds,x);
0308     h = h/sum(h);
0309     bar(x,h);
0310     xlabel('ISI');
0311     xlim([0 10]);
0312     ylabel('Count');
0313     % CV and CV2
0314     binSize = 0.1;
0315     smooth = 10;
0316     for i = 1:20,
0317         cv(i) = CV(stims0,'order',i,'binSize',binSize,'smooth',smooth,'method','fixed');
0318     end
0319     cv2 = CV(stims0,'measure','cv2');
0320     subplot(2,3,6);
0321     plot(cv,'+-');
0322     hold on;
0323     plot(cv2,'r+');
0324     xlabel('nth ISI');
0325     ylabel('CV & CV2');
0326 
0327     % Stims autocorrelograms
0328     figure;
0329     [ccg,x,y] = ShortTimeCCG(stims0,'min',1,'mode','norm','smooth',[2 0]);
0330     PlotShortTimeCCG(ccg,'x',x,'y',y);
0331     clim([0 0.01]);
0332 end
0333
FieldPSP

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
