


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 % Useful numbers 0029 f = (tref - tmin) ./ (tmax - tmin); % fraction of range that is below tref 0030 bins = linspace(0, tmax, floor(tmax*Fs)); % using tmax*Fs bins means no info is lost to rounding 0031 refractory_bin = sum((bins <= tref)); % index of refractory period 0032 0033 % Convert from times to interspike intervals 0034 isi1 = diff(sort(unit1times)); 0035 isi2 = diff(sort(unit2times)); 0036 isiT = diff(sort([unit1times; unit2times])); 0037 0038 % Calculate histograms for the intervals below tmax. 0039 isi1hist = hist(isi1(isi1 < tmax), bins); 0040 isi2hist = hist(isi2(isi2 < tmax), bins); 0041 isiThist = hist(isiT(isiT < tmax), bins); 0042 0043 % Count total # intervals below tmax (set 0 counts to 1 as a courtesy to the next step) 0044 isi1count = max(sum(isi1hist),1); 0045 isi2count = max(sum(isi2hist),1); 0046 isiTcount = max(sum(isiThist),1); 0047 0048 % Make cdfs from histograms 0049 cdfs = zeros(3, length(bins)); 0050 cdfs(1,:) = cumsum(isi1hist) ./ isi1count; 0051 cdfs(2,:) = cumsum(isi2hist) ./ isi2count; 0052 cdfs(3,:) = cumsum(isiThist) ./ isiTcount; 0053 0054 % 0055 % if ((isi1count == 1) || (isi2count == 1)), % if either original list had no violations 0056 % 0057 % 0058 % else 0059 % Compute the (scaled) difference between the initial cdfs and the combined cdfs. 0060 diffs = zeros(2, length(bins)); 0061 diffs(1,:) = sqrt((isi1count*isiTcount)./(isi1count+isiTcount)) .* (cdfs(3,:) - cdfs(1,:)); 0062 diffs(2,:) = sqrt((isi2count*isiTcount)./(isi2count+isiTcount)) .* (cdfs(3,:) - cdfs(2,:)); 0063 0064 % Max value of this difference in the region shorter than the refractory period 0065 0066 dstats(1) = max(diffs(1, 1:refractory_bin)); 0067 dstats(2) = max(diffs(2, 1:refractory_bin)); 0068 0069 % based on Kolmogorov-Smirnoff statistics (see Fee, Mitra, Kleinfeld, 1996) 0070 lambda = 0:0.01:2; 0071 pdfDstat = (0.5 * (1 + erf(lambda./(sqrt(2*f*(1-f)))))) - ... 0072 (0.5 .* exp(-2.*lambda.^2) .* (1 - erf((1-2*f).*lambda./(sqrt(2*f*(1-f)))))); 0073 cutoff = lambda(max(find(pdfDstat < 0.95))); 0074 0075 % Does either d-statistic exceed the statistical cutoff? 0076 if (any(dstats > cutoff)), allow = 0; 0077 else allow = 1; 0078 end 0079 % end 0080 0081 % The cdf at the tref bin is the fraction of spikes below tref; to get a score, 0082 % we normalize by dividing by f from above. 0083 scores = (cdfs(:, refractory_bin)) ./ f;