Biol. Cybernetics 22, 61--71 (1976) @ by Springer-Verlag 1976

Sequential Interval Histogram Analysis of Non-stationary Neuronal Spike Trains* Arthur C. Sanderson and Benjamin Kobler** Biotechnology Program and Electrical Engineering Department, Carnegie-MellonUniversity, Pittsburgh, Pennsylvania,USA Received: July 3, I975

Abstract The spike interval histogram, a commonly used tool for the analysis ofneuronal spike trains, is evaluated as a statisticalestimator of the probability density function (pdt) of interspike intervals. Using a mean square error criterion, it is concluded that a Parzen convolution estimate of the pdf is superior to the conventional histogram procedure. The Parzen estimate using a Gaussian weighting function reduces the number of intervals required to achieve a given error by a factor of 5-10. The Parzen estimation procedure has been implementedin the sequential interval histogram (SQIH) procedure for analysis of non-stationary spike trains. Segments of the spike train are defined using a moving window and the pdf for each segmentis estimated sequentially.The procedure which we have found most practical is interactive with the user and utlizes the theoretical results of the error analysis as guidelines for the evolution of an estimation strategy. The SQIH procedure appears useful both as a criterion for stationarity and as a means to characterize non-stationary activity.

purpose computers, with the prevalence of the versatile minicomputer, improved estimation techniques are preferable. A procedure for implementation of the Parzen estimation technique as an alternative to the interval histogram for the estimation of the probability density function o f interspike intervals is presented. Improved estimation of the interval pdf means, in effect, that we require fewer samples, and therefore a shorter length of data record, to obtain a given desired estimation accuracy. The sequential interval histogram (SQIH) technique is based upon this result. A shorter record length improves the temporal resolution of a density function which is changing with time. The S Q I H technique is proposed as a useful method to characterize non-stationary spike train activity.

I. Introduction The spike interval histogram (SIH) is one of the basic statistical tools used for the analysis of neuronal spike train activity. Despite the widespread use of the SIH, there has been very little attention devoted to its performance as a statistical estimator. In this paper we will reexamine the SIH with the goals of evaluating its performance and of providing guidelines for efficient estimation of the spike interval probability density function (pdf). We address this problem by calculating the expected error from three different estimation techniques including the interval histogram. It will be shown that the conventional histogram estimate is the poorest of these techniques. This result suggests that while the histogram may have had advantages in computational simplicity for special * Portions of this work were presented at the Symposium on Computer Technology in Neuroscience Research, West Virginia University Medical Center, Morgantown, West Virginia, USA, April, 1975. ** Present Address: Martin-Marietta Corporation, Baltimore, Md., USA 21227.

II. Estimation of the Interval Probability Density Function Although some applications appeared earlier, m a n y of the techniques of spike train analysis currently in use are rooted in the description of the stochastic point process formalized by Cox and Lewis (1966) (Cox and Smith, 1953; Hagiwara, 1954; Gerstein and Kiang, 1960; Perkel et al., 1967a, b; Moore et al., 1970; Gerstein and Perkel, 1972; Eckhorn and P6ppel, 1972; Bryant et al., 1973; French and Holden, 1971). In addition, several texts which treat the theory and application of stochastic point processes have appeared recently (Lewis, 1972; Srinivasan, 1973; Murthy, 1974; Cinlar, 1975). The SIH technique is commonly used in neurophysiology to characterize the temporal pattern of spontaneous nerve spike discharge. The SIH displays the relative frequency of occurrence of interspike interval lengths. A procedure for compiling a SIH is outlined in Fig. 1: (a) A "stationary" segment of a spike

62 (a)

t t~ I

t 2 t 3 t~ I II

I

t5

t6

I

I

T (d

Normalized Histogram

Z,/D t

Intervals

(b)

Spikes

t7 I

/

I

I

x = x2 = X3= x4 = X5= X6 =

(C)

where F(z) is the interspike interval distribution function,

t2 -t 1 t3 -t 2 t4-t 3 ts-tz, t6-t 5 t7 - t 6

,

F(z) = ~; f(u)du.

Histogram

3/D 1

2J~

2

o r-~,JH~-~fz/q

o

A

I

i

i

i

I

0 A 2A 3A 4A 5A 6JA

Fig. la~l. Procedure for compiling a conventional spike interval histogram. (a) A segment of a spike train of length Tcontaining N spikes (in this case, N = 7) is selected. The times of occurrence of the spikes are q , i = l , 2 ..... N. (b) Successive interval lengths, Xj, j = 1, 2..... N - 1, are determined. (c) Intervals are grouped according to their lengths. The number of intervals within a certain range, kA--+(k+I)A, is counted and plotted as a vertical bar in the histogram. (d) The total area of the histogram is normalized to one by dividing the vertical scale by D=(N-1)A. The normalized histogram is an estimate of the probability density function of interspike intervals

train record of length T seconds including N spikes is selected. (b) The N-1 interval lengths between spikes are measured. (c) The number of intervals, rig, in a given range, (kA, (k+ 1)A), k=0, 1, 2..... where A is the bin width of the histogram, is plotted versus kA, the interval length. (d). The plot is normalized (unit area) by dividing the vertical scale by (N-1)A. Implementation of the SIH procedure for a given record length, T, requires the choice of A. In practice, A is chosen largely on heuristic grounds, based on experience and an expectation of the type of histogram to be calculated. As will be shown below, guidelines can be found which aid in this choice. In statistical terms, the SIH is an estimator of the pdf of interspike intervals, f(z), where

f(z) = Prob [-event in (z, ~ + dz)] event at z = 0 and no events in between].

(3)

Although the interval histogram has been the principal tool used in neurophysiology for estimation of the pdf, considerable attention has been devoted in the statistical literature to the generalized estimation problem. In particular,..a generalized estimation, or smoothing, procedure has been discussed by Rosenblatt (1956), Whittle (1958), Parzen (1962), and Bartlett (1963). We will refer to this technique as Parzen estimation and will discuss it in detail below. A number of other techniques for estimation of both pdfs and spectral densities have been discussed (Watson and Leadbetter, 1963; Tarter and Kronmal, 1970; Good and Gaskins, 1971; Boneva et al., 1970; Fellner, 1974; Rosenblatt, 1975; Wahba, 1975a, b; Wahba and Wold, 1975a, b). The generalized form of the Parzen estimate of the pdf can be described as follows. Given N observations, Xl,X 2 XN, drawn from a population with pdff(x), the Parzen estimate, f(x), is: ....

f(x) = ~ Wx(U) dn(u),

(4)

= ~,N=I Wx(XyN,

(5)

where

N(x)-- number of observations with value less than x, and

n(x) = N(x)/N, the sample distribution function, wx(u) is a weighting function which is to be determined. First, we will show that the SIH estimate, f~i(x), which was described above is a special case of the Parzen estimate. Choose, wx(u) = 1/A, =0,

otherwise,

where

(1)

d = min (Ix k - xl). sgn (x k - x).

(8)

k

If we utilize enough observations, we would expect the relative frequency estimate to converge as follows:

sgn (x) is the signum function defined by: sgn(x) = + 1

lim

[n(N_~kl)A]=

l~(k~+l)~f(u)du

=-1

N~oo

F((k + 1)d) - F(kd) d

(2)

x > 0,

x