


ISIQUALITY Computes statistical measures of refactory period quality.
[allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs)
Returns a boolean specifying whether the refractory period statistics are
statistically worsened by combining the two input lists, along with two
outputs that give more details:
INPUTS:
unit1times, unit2times : (sec) original spike time lists
tmin : (sec) minimum possible interval
tmax : (sec) max interval for statistical comparison
tref : (sec) estimate of the refractory period
Fs : (Hz) data sampling frequency
OUTPUTS:
allow : (0 or 1) is it ok to combine the two lists?
scores : (3 x 1) isi scores for [unit1 unit2 combined] lists;
describe the normalized fraction of spikes violating
the putative refractory period (1.0 means that the
mean density of intervals below tref is the same as the
mean density from tref to tmax.
cdfs : (3 x m) isi cdfs for unit1, unit2 and combined lists,
covering the range from 0 to tmax with bin size
1/Fs. (i.e., m is floor(tmax*Fs)).

0001 function [allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs) 0002 0003 % ISIQUALITY Computes statistical measures of refactory period quality. 0004 % [allow, scores, cdfs] = isiQuality(unit1times, unit2times, tmin, tmax, tref, Fs) 0005 % 0006 % Returns a boolean specifying whether the refractory period statistics are 0007 % statistically worsened by combining the two input lists, along with two 0008 % outputs that give more details: 0009 % 0010 % INPUTS: 0011 % unit1times, unit2times : (sec) original spike time lists 0012 % tmin : (sec) minimum possible interval 0013 % tmax : (sec) max interval for statistical comparison 0014 % tref : (sec) estimate of the refractory period 0015 % Fs : (Hz) data sampling frequency 0016 % 0017 % OUTPUTS: 0018 % allow : (0 or 1) is it ok to combine the two lists? 0019 % scores : (3 x 1) isi scores for [unit1 unit2 combined] lists; 0020 % describe the normalized fraction of spikes violating 0021 % the putative refractory period (1.0 means that the 0022 % mean density of intervals below tref is the same as the 0023 % mean density from tref to tmax. 0024 % cdfs : (3 x m) isi cdfs for unit1, unit2 and combined lists, 0025 % covering the range from 0 to tmax with bin size 0026 % 1/Fs. (i.e., m is floor(tmax*Fs)). 0027 0028 % Last Modified By: sbm on Thu Jun 2 17:22:36 2005 0029 0030 % Useful numbers 0031 f = (tref - tmin) ./ (tmax - tmin); % fraction of range that is below tref 0032 bins = linspace(0, tmax, floor(tmax*Fs)); % using tmax*Fs bins means no info is lost to rounding 0033 refractory_bin = sum((bins <= tref)); % index of refractory period 0034 0035 % Convert from times to interspike intervals 0036 isi1 = diff(sort(unit1times)); 0037 isi2 = diff(sort(unit2times)); 0038 isiT = diff(sort([unit1times; unit2times])); 0039 0040 % Calculate histograms for the intervals below tmax. 0041 isi1hist = hist(isi1(isi1 < tmax), bins); 0042 isi2hist = hist(isi2(isi2 < tmax), bins); 0043 isiThist = hist(isiT(isiT < tmax), bins); 0044 0045 % Count total # intervals below tmax (set 0 counts to 1 as a courtesy to the next step) 0046 isi1count = max(sum(isi1hist),1); 0047 isi2count = max(sum(isi2hist),1); 0048 isiTcount = max(sum(isiThist),1); 0049 0050 % Make cdfs from histograms 0051 cdfs = zeros(3, length(bins)); 0052 cdfs(1,:) = cumsum(isi1hist) ./ isi1count; 0053 cdfs(2,:) = cumsum(isi2hist) ./ isi2count; 0054 cdfs(3,:) = cumsum(isiThist) ./ isiTcount; 0055 0056 % 0057 % if ((isi1count == 1) || (isi2count == 1)), % if either original list had no violations 0058 % 0059 % 0060 % else 0061 % Compute the (scaled) difference between the initial cdfs and the combined cdfs. 0062 diffs = zeros(2, length(bins)); 0063 diffs(1,:) = sqrt((isi1count*isiTcount)./(isi1count+isiTcount)) .* (cdfs(3,:) - cdfs(1,:)); 0064 diffs(2,:) = sqrt((isi2count*isiTcount)./(isi2count+isiTcount)) .* (cdfs(3,:) - cdfs(2,:)); 0065 0066 % Max value of this difference in the region shorter than the refractory period 0067 0068 dstats(1) = max(diffs(1, 1:refractory_bin)); 0069 dstats(2) = max(diffs(2, 1:refractory_bin)); 0070 0071 % based on Kolmogorov-Smirnoff statistics (see Fee, Mitra, Kleinfeld, 1996) 0072 lambda = 0:0.01:2; 0073 pdfDstat = (0.5 * (1 + erf(lambda./(sqrt(2*f*(1-f)))))) - ... 0074 (0.5 .* exp(-2.*lambda.^2) .* (1 - erf((1-2*f).*lambda./(sqrt(2*f*(1-f)))))); 0075 cutoff = lambda(max(find(pdfDstat < 0.95))); 0076 0077 % Does either d-statistic exceed the statistical cutoff? 0078 if (any(dstats > cutoff)), allow = 0; 0079 else, allow = 1; 0080 end 0081 % end 0082 0083 % The cdf at the tref bin is the fraction of spikes below tref; to get a score, 0084 % we normalize by dividing by f from above. 0085 scores = (cdfs(:, refractory_bin)) ./ f;