The rational is simply to fit a mathematical model to the
ISI distribution that includes a component to model the
refractory period.  If a neuron is well isolated, it must              
exhibit some refractory period.  If instead, it is multi-
unit, then it will exhibit instead something similar to  
a poisson process in which there is no refractory period,
or at least large probability of inter-spike intervals   
at the shortest intervals.  

The code fits two functions, the first one is a gaussian
refractory period combined with a poisson process.  Let
  X be a gaussian random variable with mean TauR and standard
    deviation of SigmaR, and let
  Y be a poisson random variable with rate R, then we
    can generate random inter-spike intervals as

 z = X + Y  where z is the random delay between spikes.
 
 The random variable z will have a distribution that is
given by the convolution of the gaussian density of X   
with the poisson density of Y, and will look much like 
a gaussian at short ISI's with a peak followed by roughly
exponential decay.  The goal of the program is to fit the
parameters TauR, SigmaR, and R in order to best describe
the shape of the actual ISI distribution.  It uses a
combination of grid search (testing many possible starting
points for the three parameters) and gradient descent to 
find the best set of parameters.  The cost function optimized
is the negative log-likelihood of the probability of the 
observed ISI distribution, say g(i), given the ISI distribution
generated by the three parameters, f(i), which is given by
     cost = sum( g(i) * log( f(i) )).

This describes how the first model is fit.  For neurons with
well defined refractory periods, we get values for TauR (the
mean refractory delay) in the range from 1.5 - 7.5 ms typically.
When we have a multi-unit, the ISI distribution is very similar
to appearing poisson (no gaussian at short ISI) and TauR is thus
best fit as being at or near 0ms.  In our dataset, and unit with
a TauR less than 1ms was double checked to confirm it was well 
isolated in the Plexon Offline Sorter.  Units with such short
TauR's indeed had questionable isolations, often a cluster that
was overlapping with another cluster, and were thus excluded
from further analysis.

If you are confident that all units in your dataset are well  
isolated from offline waveform analysis, then the ISI analysis
should not be necessary.  This ISI analysis is simply a secondary
way of validating that units included were isolated.  You
could also, just compute the ISI distributions, and look at them
by eye to confirm that in each case you see a well defined refractory
period.  It is important to compute the ISI distribution with enough
time precision to be able to see short refractroy periods on the
order of 1.5 ms, which is physiologically possible, especially among
neurons with bursty firing properties.  So if you discretize the
ISI distribution in bins of about 0.5ms this would probably be
good enough.  I think the AV Herz paper discretized at 0.1ms. 

Lastly, the second model that is fit is able to deal with neurons
that have bursty firing properties (in which there is a peak at
very short ISI's corresponding to the intervals between spikes
during a burst events).  This model may provide a better fit when
neurons have significant bursting.  Part of the reason we have
developed these models is not really for validating unit isolations,
but more specifically for fittin ISI distributions and assessing
how aspects of the ISI distribution, such as burstiness, change
with attention.

Question 1):  are you assuming the firing rate is constant over
the interval in which the ISI distribution is analyzed?

That is absolutely correct.  In our task a stimulus gradually
enters into the visual receptive field and remains there for
a 1000ms period.  We only use the last 800ms of this period,  
which is relatively free of response onsets, for the ISI analysis.
Even so, the assumption that rate is stable in this 800ms is
probably not completely correct.  It appears stable when you
compute the PSTH over trials, but in other analysis, we can
clearly see that firing rate is modulated by local field osciallations
during this 800ms period trial by trial, and thus a more complete
model will ultimately be necessary that captures fluctations in rate
due to osciallations in the local field as well.

You may still be able to fit ISI distributions over the entire
trial.  This would probably still give you a good estimate of
the refractory period, although the rate parameter of the poisson
would not be as meaningful anymore.

Question 2):  what is the format of the spike data structures provided?

You are correct, each row represents at trial with 1600 values,
and this corresponds to 800ms with time discretization of 0.5ms.
Each value will be 0 (no spike) or 1 (spike).  There should be
roughly 70-100 rows, each corresponding to one trial.


   Jude