## Abstract

The ability of an animal to detect weak sensory signals is limited, in part, by statistical fluctuations in the spike activity of sensory afferent nerve fibers. In weakly electric fish, probability coding (P-type) electrosensory afferents encode amplitude modulations of the fish's self-generated electric field and provide information necessary for electrolocation. This study characterizes the statistical properties of baseline spike activity in P-type afferents of the brown ghost knifefish, *Apteronotus leptorhynchus.* Short-term variability, as measured by the interspike interval (ISI) distribution, is moderately high with a mean ISI coefficient of variation of 44%. Analysis of spike train variability on longer time scales, however, reveals a remarkable degree of regularity. The regularizing effect is maximal for time scales on the order of a few hundred milliseconds, which matches functionally relevant time scales for natural behaviors such as prey detection. Using high-order interval analysis, count analysis, and Markov-order analysis we demonstrate that the observed regularization is associated with memory effects in the ISI sequence which arise from an underlying nonrenewal process. In most cases, a Markov process of at least fourth-order was required to adequately describe the dependencies. Using an ideal observer paradigm, we illustrate how regularization of the spike train can significantly improve detection performance for weak signals. This study emphasizes the importance of characterizing spike train variability on multiple time scales, particularly when considering limits on the detectability of weak sensory signals.

- electrosensory afferent
- electrolocation
- interspike interval analysis
- Markov process
- spike train variability
- weak signal detection

Survival in an animal's natural environment is dependent on the ability to detect behaviorally relevant stimuli, such as those caused by predators and prey. Being able to reliably and efficiently detect such signals at weak levels confers a competitive advantage. Thus many sensory systems, including the electrosensory system discussed here, have presumably experienced selective pressures over the course of evolution to improve detection performance for weak sensory signals.

The decision of whether or not a stimulus is present must ultimately be based on a change in the spike activity of primary afferent nerve fibers. In many cases, this change must be detected in the presence of ongoing spontaneous activity. Intuitively, a subtle change in spike activity caused by a weak external signal should be easier to detect when the baseline activity is regular and predictable than when it is irregular and subject to random fluctuations. To understand the limits on signal detection performance, it is thus important to characterize the variability of baseline activity in primary afferent spike trains.

A common approach for characterizing spike train variability is by analysis of the first-order interspike interval (ISI) distribution (Hagiwara, 1954; Moore et al., 1966;Ratliff et al., 1968; Gabbiani and Koch, 1998). The coefficient of variation (the SD divided by the mean) of the ISI distribution provides a convenient measure of the variability in the arrival time between successive spikes. It is important to realize, however, that the first-order ISI distribution only provides information about short-term variability on time scales comparable to the mean ISI. Long-term variability, over time periods containing multiple spikes, must be measured using other techniques, such as analysis of higher-order interval distributions (Rodieck et al., 1962; Moore et al., 1966) or spike count distributions (Barlow and Levick, 1969a,b; Teich and Khanna, 1985).

Measurements of long-term variability for primary afferent spike trains are not commonly encountered in the literature. These measures would indeed be redundant if afferent spike activity could be adequately modeled as a renewal process. For a renewal process, successive intervals in the ISI sequence are independent and identically distributed (Cox, 1962) and therefore, higher-order interval and count distributions can be derived knowing only the first-order ISI distribution. Thus, variability on all time scales can be computed. However, when spike activity arises from a nonrenewal process there will be correlations and history-dependent effects in the ISI sequence. In such cases, the first-order ISI distribution does not provide sufficient information to predict long-term spike train variability or to set limits on signal detection performance.

In this paper, we analyze the variability of baseline spike activity recorded from P-type (probability coding) electrosensory afferent fibers in the weakly electric fish *Apteronotus leptorhynchus*(brown ghost knife fish). Objects near the fish that differ in impedance from the surrounding water modulate the self-generated electric field because of the fish's electric organ discharge (EOD). These modulations provide sensory cues that allow the fish to hunt and navigate in the dark using electrolocation (Rasnow, 1996) (for review, see Bullock and Heiligenberg, 1986). P-type afferents respond to the strength of amplitude modulations (AMs) by increasing or decreasing their probability of firing (Scheich et al., 1973; Bastian, 1981; for review, see Zakon, 1986). Their AM response characteristics have been well studied (Hagiwara et al., 1965; Scheich et al., 1973; Hopkins, 1976; Bastian, 1981; Shumway, 1989; Wessel et al., 1996; Xu et al., 1996; Nelson et al., 1997), but variability of baseline spike activity has not been fully characterized.

Infrared video recordings of prey capture behavior in*Apteronotus* performed in our laboratory have been used to estimate the behavioral threshold for detecting small prey (*Daphnia magna*, 2–3 mm in length) in the dark. At the time of detection, we estimate that the prey gives rise to an AM signal that transiently changes the firing probability of P-type afferents by only ∼1% (Nelson and MacIver, 1999). In this study we show that P-type afferent spike trains are irregular as judged by the first-order ISI distribution but that there is an underlying nonrenewal process that serves to make the spike train more regular over longer time intervals. In an ideal observer framework (Green and Swets, 1966) this regularity effectively reduces the detection threshold for weak stimuli, such as those caused by small prey.

## MATERIALS AND METHODS

#### Electrophysiology

Extracellular recordings were made from isolated P-type afferent fibers of weakly electric fish *Apteronotus leptorhynchus.*Surgical and nerve recording procedures used here are identical to those described in an earlier study (Xu et al., 1996). Briefly, fish were anesthetized in 100 ppm tricaine methanesulfonate (MS-222; Sigma, St. Louis, MO) and immobilized with an intramuscular injection of 3 μl 10% gallamine triethiodide (Flaxedil; Sigma). P-type afferent fiber activity was recorded from the posterior branch of the left anterior lateral line nerve (pALLN), which innervates trunk electroreceptors. Recordings were made in the presence of the EOD of the fish with no other stimulus present. We refer to activity under these conditions as “baseline” activity, in contrast to spontaneous activity that would be obtained if the EOD were silenced. Action potentials were recorded from individual pALLN fibers with glass microelectrodes (impedance of 10–30 MΩ) filled with 3 mKCl solution. Neural activity and EOD waveforms were digitized at 17 kHz and stored for offline analysis. Spike events in the nerve recording were identified by a threshold criteria and time-stamped with a resolution of 60 μsec. All data analysis was performed on Sun workstations using custom software and the MATLAB programming environment (The MathWorks).

#### Spike train representation

*Apteronotus leptorhynchus* has a continuous quasi-sinusoidal EOD waveform with a fundamental frequency f that ranges from 750 to 1000 Hz depending on the individual. P-type units fire at most once per EOD cycle and randomly skip cycles between successive spikes. On average, a typical unit fires on about one-third of the EOD cycles. This ratio is referred to as the per-cycle probability of firing p. In the presence of a stimulus, the per-cycle probability is modulated by stimulus intensity, and hence P units are called probability coders.

When spike times are sampled at intervals smaller than the EOD period, the timing of spikes within the cycle can be observed. This is illustrated in the ISI histogram shown in Figure1*A.* The peaks of the ISI distribution are separated by one EOD period. The width of each peak reflects the variability in firing phase within the EOD cycle. In subsequent analyses in this paper, we only keep information about the occurrence of a spike in an EOD cycle and discard information about the phase within the cycle. That is, we resample the spike train at the EOD frequency f. The effect of this resampling is seen in Figure1*C*, which is the discrete time ISI histogram for the fiber shown in Figure 1*A.* Time is measured in EOD cycles and hence, ISIs assume only integer values j ≥ 1, where j is the number of cycles to the next spike. For the remainder of this work, we only consider discrete-time spike trains sampled at the EOD rate.

The resampled spike train is a realization of a discrete time stochastic process x(i) where i ≥ 1 is the EOD cycle number. Furthermore, x(i) = 1 if there is a spike in cycle i and is zero otherwise. The mean per-cycle firing probability is given by p = n/T where n = ∑_{i=1}^{T} x(i) is the total number of spikes observed in a record of duration T EOD cycles. The stochastic process x(i) can be characterized in terms of the statistical properties of the time intervals between spikes (interval analysis) or by the statistical properties of spike counts in time windows of fixed durations (count analysis) (Cox and Lewis, 1966).

#### Interval sequences and distributions of order k

Let t_{i} represent the EOD period number in which the i th spike occurs. We set the time origin to be such that t_{1} = 1. Interval sequences of order k are defined as (Rodieck et al., 1962;Moore et al., 1966)
Equation 1where n is the total number of spikes. The sequence S_{k} is a discrete-time stochastic process of strictly positive integers. The first-order interval sequence S_{1}(i) = t_{i+1} − t_{i} is the sequence of ISIs, S_{2}(i) is the sequence of times between every second spike, etc. Let I_{k} be the normalized k th order interval histogram (IH) constructed from S_{k}. Then I_{k}(j) is the probability of observing an interval of length j in the sequence S_{k}. We denote the mean and variance of I_{k} by
_{k} and Var (I_{k}), respectively. In this paper, interval sequences and IHs were calculated for orders k up to 4096. In the neurophysiology literature, I_{1}is used extensively and is referred to as the ISI distribution. The shape and statistical properties (mean and variance) of the ISI distribution are often used to characterize firing patterns and variability of spike timing. However, the ISI provides a characterization of spike variability only on time scales comparable to the mean ISI. Higher-order interval distributions provide information about variability over longer time scales and about dependencies in the ISI sequence.

Dependency of an interval on the immediately preceding interval was analyzed using the joint interval histogram I(j_{1}, j_{2}), which reflects the probability of observing an ISI of length j_{1} followed by an ISI of length j_{2} (Rodieck et al., 1962). The joint interval histogram is constructed from the sequence S_{1} by binning all overlapping tuples (S_{1}(i), S_{1}(i + 1)), where 1 ≤ i ≤ n − 1. Figure 1, *B* and*D*, shows plots of the joint interval histogram corresponding to the sampling resolutions of Figures 1*A* and *C*, respectively.

#### Count distributions

An alternate analysis of spike train variability can be made using spike count distributions (Barlow and Levick, 1969a,b;Teich and Khanna, 1985). Proceeding as for interval analysis, count statistics were obtained from two measures: count sequences and count histograms. Let T be the number of EOD periods in which a count is to be performed. Then we define the count sequence N_{T} by counting the number of spikes occurring in blocks of T contiguous EOD cycles. This can be expressed as:
Equation 2where i refers to the block number in the sequence and n is the total number of recorded spikes. The normalized count histogram C_{T} was also calculated from N_{T}, so that C_{T}(m) is the probability of obtaining m spikes in a count window of duration T EOD cycles. The mean and variance of C_{T} are denoted by
_{T} and Var (C_{T}), respectively. Count sequences and histograms were calculated for T values ranging from 20 to as many as 50,000 EOD cycles, depending on data availability and subject to a minimum of 10 counting blocks.

#### Spike train models

To gain an understanding of the process that generates the observed spike train data x(i), we created surrogate spike trains using three model systems that reproduce certain statistical features of the afferent data.

*Binomially generated spike train (B).* The classic description of a P-type afferent is that it is a “probability coder.” Namely, the afferent fires irregularly with a per-cycle probability p that is constant under baseline conditions, but is subject to modulation in the presence of external stimuli. Thus, the simplest model of P-unit baseline activity is a process that emits a spike in any given EOD cycle with constant probability p independent of previous history. In a continuous-time framework, a constant firing probability per unit time gives rise to a homogeneous Poisson process, which frequently serves as a basis for models of spontaneous activity in neural spike trains (Tuckwell, 1988). In the discrete-time framework considered here, a constant firing probability per time step (one EOD cycle) gives rise to a binomial process. Whereas the Poisson process has exponentially distributed ISIs and Poisson-distributed spike counts, the binomial process has geometrically distributed ISIs and binomially distributed counts (see Eqs. 10 and 12, Appendix A). Surrogate binomial spike trains B were generated by shuffling the observed spike sequence x(i) to remove all dependencies between adjacent EOD periods. Per-cycle firing probability p remains unchanged because shuffling preserves the total number of spikes and EOD periods.

*Zeroth-order Markov process (M _{0}).* The binomially generated spike train matches p but does not guarantee that the interval distributions will be the same as the data. A surrogate spike train that preserves p as well as the ISI distribution I

_{1}, can be generated by randomly shuffling the ISI sequence S

_{1}, rather than shuffling the spike train x. This preserves the total number of intervals and the distribution of intervals but removes dependencies between neighboring intervals in the sequence (Longtin and Racicot, 1997a). For reasons discussed below, this process will be referred to as M

_{0}, indicating that it is a zeroth-order Markov process.

*First-order Markov process (M _{1}).* Interval distributions of experimentally observed spike trains often exhibit dependencies on prior activity and are thus nonrenewal (Kuffler et al., 1957; Werner and Mountcastle, 1963;Teich et al., 1990; Lowen and Teich, 1992). Surrogate spike trains that preserved dependencies between adjacent ISIs were constructed from the afferent data. First, all adjacent pairs of intervals (j

_{i}, j

_{i+1}) were tabulated and sorted into groups, each having identical first element. A group, say with first element j

_{a}, was selected at random, and a tuple was drawn, say (j

_{a}, j

_{b}). The next tuple was drawn at random from the group which had first element j

_{b}. If this tuple was (j

_{b}, j

_{c}), the resultant of the two draws was the triplet (j

_{a}, j

_{b}, j

_{c}). All draws were made without replacement. Continuing in this manner, an ISI sequence was constructed that had joint probability I(j

_{1}, j

_{2}) that matched the data. As discussed below, the resulting sequence is a first-order Markov process and will be denoted by M

_{1}.

The binomial (B) and zeroth-order Markov (M_{0}) processes are examples of renewal processes (Cox, 1962) because successive intervals in the ISI sequence S_{1} are independent and identically distributed. The first-order Markov process (M_{1}) is a nonrenewal process because intervals are not independent. Appendix summarizes some results for renewal processes.

#### Measures of variability

A common measure of spike train variability is the coefficient of variation of the ISI distribution CV_{I}, defined as the SD of the ISI distribution divided by its mean. Because it is a dimensionless quantity, it can be used for comparing the variability of two distributions even when they differ in their means. To measure variability of ISIs on different time scales, we also computed CV_{I} for higher order interval distributions. The coefficient of variation CV_{I}(k) for the k th order interval distribution I_{k}is:
Equation 3

Another useful measure is the variance-to-mean ratio of I_{k}, denoted F_{I}(k):
Equation 4

Although this measure has the disadvantage that it is not dimensionless, it has the useful property that for all orders k, it is constant for a renewal process (see Appendix ). For any process that is more regular than a renewal process, the variance-to-mean ratio decreases with increasing k.

Analogous to the measures of variability for interval sequences, it is possible to define similar measures for count sequences N_{T}. Proceeding as above, the coefficient of variation CV_{C}(T) for the count distribution C_{T}(i) is defined as:
Equation 5

For spike count distributions, the variance-to-mean ratio is called the Fano factor (Fano, 1947). It is denoted F_{C}(T) and defined as:
Equation 6

#### Correlation analysis

History-dependent effects were analyzed by considering the serial correlation coefficient (SCC) ρ_{l} of the first-order sequence S_{1}, where l is the lag in terms of the number of intervening intervals. The ρ_{l} were computed from:
Equation 7where j_{1}, … , j_{M}is a sequence of M consecutive ISIs that can start anywhere in the ISI sequence S_{1}. The values of ρ_{l} range from +1 (perfect correlation) to −1 (perfect anti-correlation), and ρ_{l} = 0 when intervals are uncorrelated.

To determine if the ρ_{l} were significantly different from zero, the sequence S_{1} was divided into nonoverlapping blocks each having M = 1000 elements S_{1}(i), … , S_{1}(i + M − 1), and the ρ_{l} were determined for each block. The sequence S_{1} was then shuffled to eliminate dependencies between intervals, if any, and the ρ_{l} were again evaluated for the same number of blocks M. For each lag, the unshuffled and shuffled SCCs were tested under the null hypothesis that the two populations were identical. A Wilcoxon rank sum test was used to test the hypothesis at a significance level of *p* = 0.01.

#### Analysis of Markov order

SCCs do not completely characterize dependencies between ISIs. This can be seen from Equation 7 where the SCC at lag l depends only on the pairs (S_{1}(i), S_{1}(i + l)) but does not depend on the intervening (l − 1) lags. Higher-order history-dependent effects in the interval sequence S_{1} can be modeled by Markov chains (Nakahama et al., 1972; van der Heyden et al., 1998). The ISI sequence S_{1} can be described by a Markov chain of order n, if intervals in the chain are dependent on exactly n preceding intervals. If the intervals are independent, the process is referred to as zero-order Markov (M_{0}). Examples of such processes are the binomially generated spike train (B) and shuffled ISI process M_{0} described above. The process M_{1}, on the other hand, is first-order Markov because the statistics of the current interval are completely determined given the previous interval.

Consider a sequence of ISIs {j_{m}, j_{m−1}, … , j_{0}}, where j_{r} refers to the r th lag relative to the current ISI j_{0}. For a Markov chain we can define the m th order transition probability as p(j_{0}‖j_{m}, … , j_{1}), which is the probability of observing the interval j_{0} given that we have observed the sequence {j_{m}, j_{m−1}, … , j_{1}} in the immediate past. Note that the definition of the m th order transition probability does not say anything about the order of the chain itself. The Markov order n of the chain is defined as the smallest value of n for which p(j_{0}‖j_{m}, … , j_{n}, … , j_{1}) = p(j_{0}‖j_{n}, … , j_{1}) for all m ≥ n. Hereafter, we use the symbol m to denote the order of the transition probability, and the symbol n as the fixed number representing the order of the Markov chain. The transition probabilities p(j_{0}‖j_{m}, … , j_{1}) can be estimated from the experimentally observed ISI sequence by counting all occurrences of j_{0} immediately after the tuple (j_{m}, … , j_{1}).

Given an experimentally observed ISI sequence, we wish to determine the Markov order of the underlying process that generated the sequence. This can be done by comparing transition probabilities obtained from the data with transition probabilities of surrogate spike trains that are constructed to be of known Markov order. Statistical comparisons can be made using the m th order conditional entropy h_{m} as a test statistic (van der Heyden et al., 1998). The h_{m} are defined as:
Equation 8
Equation 9for all possible tuples (j_{m}, … , j_{0}) having nonzero joint probability p(j_{m}, … , j_{0}). The conditional entropy satisfy h_{m+1} ≤ h_{m}, for all m (Shannon, 1948). For an n th order Markov chain, h_{m} = h_{n} for all m ≥ n.

We follow the procedure of van der Heyden et al. (1998)for testing the order of a Markov chain. Hypothesis testing was performed for increasing orders m, beginning with m = 0. The null hypothesis was that the process is Markov order m. The alternative hypothesis was that the process is order (m + 1) or greater. The test statistic was the (m + 1) th order conditional entropy h_{m+1} given by Equation 9. The statistic h_{m+1} was evaluated for both afferent and surrogate data sets. A total of R_{s} = 49 surrogate data sets were generated, and the rank r of the afferent data set was determined. Because the test is one-sided, the*p* value for the test is p = r/(R_{s} + 1). The null hypothesis was rejected if*p* ≤ 0.05.

Hypothesis testing was performed for increasing orders m until the null hypothesis could not be rejected or until the testing was terminated because of insufficient numbers of surrogate sequences. The criteria for terminating the test was N_{m} < N/R_{s}, where N was the number of intervals in S_{1}, and N_{m} was the number of distinct (m + 1) tuples extracted from the data set. In this case only a lower-bound for the order n could be estimated.

#### Weak signal detection

The impact of spike train regularity on signal detection performance was estimated for afferents and the three matched spike train models. The detection algorithm was based on an ideal observer paradigm (Green and Swets, 1966) using spike count distributions. For each binary spike train x(i) (see above), signal windows of duration 100 EOD periods were selected at regular intervals of 300 EOD periods. A random offset (uniform between 0 and 99) was added to the starting position of each signal window. A trial consisted of the addition of a constant number of spikes (n_{s}) to each signal window. Spikes were randomly distributed and added only in those EOD periods that did not already contain a spike. A sequence of spike counts was then generated using Equation 2 with T = 100 EOD periods. If the count exceeded a fixed threshold (see below) a “hit” was generated for that counting window. Because the threshold can be exceeded because of random fluctuations in the baseline even when there is no signal present, a percentage of the hits will be false alarms. Let N_{s} denote the total number of counting windows where signal + baseline is present, with N_{hs} of these windows receiving a hit. Similarly, let N_{b} denote the total number of counting windows where only baseline is present, with N_{hb}of these windows receiving a hit. Then detection probability is given by P_{d} = N_{hs}/N_{s}, and false alarm probability is given by P_{fa} = N_{hb}/N_{b}. The threshold was chosen such that the false alarm probability P_{fa} was 0.001 or less, motivated by the low rate of false strikes observed in our prey capture studies (Nelson and MacIver, 1999). Trials were repeated for n_{s} = 1, 2, … , 30, and P_{d} was evaluated as a function of n_{s}.

## RESULTS

Baseline spiking activity in the absence of any external stimulus other than the fish's ongoing EOD was recorded from 52 individual P-type afferent fibers in eight fish. The EOD frequency of an individual fish was constant and ranged from 750 to 1000 Hz. Afferent baseline record lengths ranged from 83 to 2048 sec, with a median of 428 sec. The firing rate for individual afferent fibers was nearly constant over the duration of the recording. Baseline firing ranged from 65 to 575 spikes/sec, with a population mean of 260 ± 124 spikes/sec (mean ± SD). These values are in agreement with previously reported results for P-type afferents in this species (Xu et al., 1996).

### Interspike interval analysis

The discrete time interspike interval histogram for a representative spike train is shown in Figure 1*C.* For this unit, ISIs range from 1 to 7 EOD cycles (
_{1} = 2.9 EOD periods; CV_{I}(1) = 0.47). Mean ISIs of afferents ranged from 1.4 to 14.2 EOD cycles with a population mean of 4.2 ± 2.3 EOD cycles (Fig. 2*A*). The coefficient of variation of the ISI distribution ranged from 0.15 to 0.79 with a population mean of 0.44 ± 0.16 (Fig.2*B*). Thus, on average, the SD of the ISI is ∼44% of the mean, which reflects a considerable degree of variability in the baseline ISI distribution.

A representative joint interval histogram, which provides information about dependencies between adjacent intervals is shown in Figure1*D.* Long intervals were more likely to be followed by short intervals and vice versa. Almost all fibers demonstrated a similar pattern.

The first-order ISI provides a characterization of spike time variability only on time scales comparable to the mean interspike interval. To characterize spike variability over longer time scales, the statistical properties of the higher-order interval distributions I_{k} must be considered. In particular it is informative to examine how the SD (ς_{I}(k)), coefficient of variation CV_{I}(k), and the variance-to-mean ratio F_{I}(k) vary as a function of interval order k. These are shown in Figures3*A–C* respectively, for a representative afferent fiber (○, *solid line*). The ς_{I} measured in EOD periods (Fig.3*A*) grows slowly as k increases from 1 to ∼60. Thereafter it increases rapidly (note logarithmic coordinates). Given that the mean of I_{k} obeys
_{k} = k
_{1}, i.e., it grows linearly with k, the initial slow growth in ς_{I} results in a decrease in the coefficient of variation CV_{I}(k) (Fig. 3*B*) from an initial value of 0.58 (k = 1) until it reaches a plateau level of ∼0.007 for interval orders k > 100. Thus, for high-order intervals, the SD is ∼0.7% of the mean. This represents a reduction in CV_{I} on long-time scales by a factor of ∼80 relative to the CV_{I} for the ISI. Similarly, for this unit, the variance-to-mean ratio F_{I}(k) drops rapidly with increasing k from an initial value of 0.69 and reaches a minimum value of 0.02 for k = 59. Thereafter, F_{I} increases steadily for larger k. Both the minima in the F_{I} curve and the knee in the CV_{I} curve are a consequence of the transition in the behavior of ς_{I} from slowly varying to rapidly increasing (Fig. 3*A*).

The trends observed in ς_{I}, CV_{I} and F_{I} for this representative unit were consistent across the entire population. Figure 3*D–F*illustrates the population distribution for all 52 fibers by plotting the median value (○, *thick line*) bracketed by the upper and lower quartiles (*thin lines*). For each fiber, we also determined the interval order k_{min}, which gave rise to the minimum variance-to-mean ratio F_{I} (Fig.3*C*). The k_{min} were distributed across the population as shown in the gray histogram (Fig.3*F*). The k_{min} for the population had a mean value of 42 ± 35. When these values are converted to EOD periods (
_{1}k_{min}) or milliseconds (
_{1}k_{min} 10^{−3}/f), they correspond to a mean time of 152 ± 121 EOD periods (or 176 ± 141 msec). The population means for the various measures at k_{min} were: ς_{I}, 2.9 ± 1.4 EOD periods; CV_{I}, 0.03 ± 0.02; and F_{I}, 0.08 ± 0.07 EOD periods. For 24 of 52 fibers, the SD was <2% of the mean interval length at k_{min} (Fig.4*A*). This should be contrasted with the SD of the ISI where the mean CV_{I}(1) for the population was 44% of the mean ISI (Fig.2*B*).

The higher-order interval measures suggest that there may be two different regimes in time demarcated by k_{min}. When k < k_{min}, ς_{I}grows very slowly, and hence, both CV_{I} and F_{I}decrease as k^{−1}. That is, the interval distribution becomes less variable as k increases. When k > k_{min}, ς_{I}grows in proportion to k, and this causes CV_{I} to plateau and F_{I} to increase as k. In this regime, variability increases with increasing k. To analyze and interpret these results, interval statistics for afferents were compared with surrogate spike trains B, M_{0}, and M_{1} described earlier.

### Comparison with renewal process models

#### Binomial model

For each afferent fiber, a binomial model B was constructed from the p (per-cycle probability of firing) value of the fiber. Under baseline conditions p is a constant and is independent of firing activity in preceding EOD cycles. This model follows from a description of P units as probability coders (see Materials and Methods). The mean p for the population of fibers was 0.31 ± 0.14 and agrees with previously published results (Xu et al., 1996).

ISI and joint interval distributions for a representative fiber are shown in Figure 5, *A1* and*A2* (same data as in Fig. 1*C,D*). The corresponding distributions for a binomial model with the same p value are shown in Figure 5, *B1* and *B2.* Although the binomial model has the same *p* value, and thus the same mean ISI
_{1} = 1/p, it does not match the ISI or joint interval distributions of the data.

Nevertheless, the binomial model serves as a useful benchmark because the coefficient of variation CV_{I} and variance-to-mean ratio F_{I} can be computed analytically (see Appendix ). For the k th-order interval distribution, ς_{I}(k) =
/p, CV_{I}(k) =
, and F_{I}(k) = (1 − p)/p = constant. Spike trains generated from this model were subject to the same analysis as the data. These are shown in Figure 3*A–C*(□, *dotted line*). Because ς_{I}(k) ∝ k^{1/2}, in logarithmic coordinates, ς_{I}(k) grows linearly with interval order k with a slope of 1/2 (Fig. 3*A*). It can be seen that the afferent has a much slower increase in ς_{I} for k < 60, but a much faster increase for orders >60. CV_{I}(k) for the afferent data decreases more rapidly with increasing k (nearly as k^{−1}) than the binomial model for which CV_{I}(k) ∝ k^{−1/2}(Fig. 3*B*, □) but eventually reaches a plateau for large interval orders. The distinction between the binomial model and the afferent data is most strikingly illustrated by the plot of the variance-to-mean ratio F_{I}(k) shown in Figure3*C.* For the binomial model F_{I}(k) = constant, but the afferent data shows a strong dependence on interval order, initially dropping rapidly (nearly as k^{−1}) and reaching a minimum value for an interval order near 60. The population medians of ς_{I}, CV_{I}, and F_{I} for a binomial process with mean p obtained by averaging over all afferents are also shown in Figure3*D–F* for comparison.

For each afferent fiber, we also determined the ratio of the F_{I}(k_{min}) for surrogate divided by afferent. Because the means of I_{k} are the same for afferent and the surrogate data sets (by construction), the above ratio is also the ratio of the variances of afferent and binomial model at k = k_{min}. This ratio provides a measure of how much more regular the afferents are in comparison with the binomial process. Figure6*A* shows the distribution of this number for the population. At k = k_{min} the afferents exhibited variance-to-mean ratios that were on the average 50 times smaller (mean, 50 ± 27) than the F_{I} for binomial spike trains.

#### Zeroth-order Markov model (shuffled ISI)

The binomial model agrees with the data only in the mean value of the ISI distribution, but it fails to accurately describe the ISI distribution. The zeroth-order Markov model (M_{0}) was constructed for each afferent by shuffling the ISI sequence S_{1} for that unit (see Materials and Methods). Thus, it matches the ISI distribution of the afferents. Figure 5, *C1* and *C2*, shows the ISI and joint interval distributions for this model. It can be seen that the model matches the ISI of the afferent (Fig. 5, compare *C1*,*A1*), however it does not correctly model the joint interval distribution (Fig. 5, compare *C2*, *A2*).

Because M_{0} is a renewal process, the coefficient of variation CV_{I} and variance-to-mean ratio F_{I}can be computed analytically (see Appendix ): CV_{I}(k) = CV_{I}(1)/
, where CV_{I}(k) = Var (I_{1})^{1/2}/
_{1}, and F_{I}(k) = Var (I_{1})/
_{1} = constant. Results for this model matched to a representative afferent are shown in Figure 3*A–C* (◊, *dash–dot line*). Note that ς_{I}(1), CV_{I}(1), and F_{I}(1) (i.e., the y intercepts) are the same for M_{0} and afferent because M_{0} was constructed to match the first-order statistics of the data. The population ς_{I}, CV_{I}, and F_{I} curves for this model are shown in Figures3*D–F*, respectively (*dash–dot line*).

Although M_{0} is more regular than the binomial model B, it still does not capture the intermediate-term regularity of the afferent data. At the minimum of the variance-to-mean ratio curve, F_{I} for the data was on average 14 times smaller (mean, 13.8 ± 10.1) than for M_{0}(Fig. 6*B*). Whereas this is an improvement over the binomial model, the reduction in variability with increasing order of intervals is not as great as that demonstrated by afferents.

The above results show that variability in afferent interval distributions cannot be accounted for by an underlying renewal process. In particular, at the minima in the F_{I} curve, which occurs on a time scale of ∼60 intervals (∼200 msec), afferents have a variance that is an order of magnitude smaller than that predicted by a renewal process. Thus, we must look at nonrenewal properties for a better understanding of this dramatic reduction in spike train variability.

### Comparison with nonrenewal models

Renewal models are “memoryless”, that is interval distributions are independent of previous intervals. However, this is not the case for P-type afferent data, as illustrated by the joint interval distribution (Fig. 5*A2*). Long intervals in the sequence are more likely to be followed by short intervals and vice versa. In contrast, neither of the renewal models (B, M_{0}) discussed above incorporate this dependency, as can be seen from Figure 5, *B2* and*C2.*

The simplest nonrenewal model that incorporates the observed long–short dependency is the first-order Markov process M_{1} described earlier (see Materials and Methods) where the probability of generating the next interval depends only on the current interval. Figure 5*D2* shows the joint interval histogram for the fiber shown in Figure 5*A2.* By construction, the ISI distribution I_{1}(j_{1}) (Fig. 5*D1*), joint interval distribution I(j_{1}, j_{2}) (Fig. 5*D2*), and second-order interval I_{2}(j_{1}) distribution (data not shown) are identical to those of the afferent.

Results for the M_{1} model are shown in Figure3*A–C* (▵, *dashed line*). ς_{I}, CV_{I}, and F_{I} are identical for M_{1} and afferent for k = 1 and k = 2 (by construction), but thereafter, M_{1} does not fully capture the continued decrease in variability observed in the afferent data, although it is better than the two renewal models B and M_{0}.

At the minimum order k_{min}, F_{I} for the afferents was on average a factor of four times smaller than the M_{1} model (mean, 4 ± 2) (Fig.6*C*). Thus, even a nonrenewal model that correctly incorporates the long–short correlations exhibited in the data cannot adequately account for the longer term regularity of afferent spike trains. This result indicates that dependencies in the underlying nonrenewal process go beyond the immediately preceding ISI. The analyses in the following sections are intended to address how far back in time the “memory” effect extends.

### Interval correlations

Correlations over longer durations were analyzed for each afferent and its matched M_{1} model using serial correlation coefficients (Eq. 7). Figure7*A* shows the population distribution of the correlation coefficient between adjacent ISIs (ρ_{1}). The ρ_{1} ranged from −0.23 to −0.82 with a population mean of −0.52 ± 0.14. That is, all afferents exhibited strong long–short ISI dependency. The distribution of correlation coefficients in the population for lags of 2 and 3 are shown in Figure 7, *B* and *C.* Most afferents exhibit positive correlations for lag = 2 (
_{2}, 0.10 ± 0.18) and weak negative correlations for lag = 3 (
_{3}, −0.07 ± 0.11). However, for lags larger than one, the coefficients become smaller, and we used more rigorous statistical criteria to determine the extent and significance of the ρ_{k}.

Serial correlation coefficients were estimated for all fibers for the first 10 lags. Coefficients are shown for three representative fibers in Figure 8 (*solid line*), and the first-order Markov model M_{1} (*dash-dot line*). Spike trains from the renewal models B and M_{0} did not exhibit any correlations (i.e., ρ_{k} = 0, for all k) because of independence of ISIs (data not shown). For both afferent and model M_{1} we used a shuffled data technique to determine the maximum lag k for which the ρ_{k} were significantly different from zero (see Materials and Methods).

Although all fibers exhibited statistically significant correlations extending back several lags, so did the M_{1}model. In fact, in almost all cases, the M_{1}model showed much stronger correlations than the afferent for all lags k ≥ 2. Thus, even a Markov process that is first-order (M_{1}) can have nonzero SCCs that extend much further back in time than a single interval.

The following example illustrates the converse effect that the SCC for lag l can be small even when there are strong l th order dependencies in the data. To illustrate this point we took the ISI sequence S_{1} for the representative afferent illustrated in Figure 1 and analyzed the distribution of Y = S_{1}(i), given X = S_{1}(i − 4) for two example patterns: (X, ?, ?, ?, Y) and (X, 2, 2, 2, Y). In the first pattern the “?” represents any arbitrary ISI, i.e., X and Y were not conditioned on the intervening intervals (S_{1}(i − 3), S_{1}(i − 2), S_{1}(i − 1)). However, in the second sequence, X and Y were conditioned on the occurrence of a specific sequence (S_{1}(i − 3) = 2, S_{1}(i − 2) = 2, S_{1}(i − 1) = 2). We then calculated the joint probability of Y given X as for Figure 1*D* (see Materials and Methods). The results are shown in Figure 9.

Figure 9*A* shows that X and Y are only weakly correlated (ρ_{4} = 0.014; *p* = 0.01). It is useful to compare this histogram with that of the zeroth-order Markov process for the same afferent as shown in Figure5*C2.* The zeroth-order model is a renewal process where intervals are uncorrelated for all lags k. The similarity between the two histograms suggests a loss of correlation. In contrast, the pattern (X, 2, 2, 2, Y) has a joint distribution shown in Figure 9*B*, and it can be seen that conditioning very strongly influences the correlation structure of the joint distribution. Because serial correlation analysis is based on unconditioned distributions, long-term memory effects may average out.

### Analysis of Markov order

The above examples illustrate the shortcomings of serial correlation analysis of memory effects and indicate that a more direct test of Markov order is needed. Toward this end we applied a statistical test to explicitly test the Markov order of the process (see Materials and Methods). Each afferent was tested against surrogate data sets that were constrained to match the afferent up to a specified order m. The test statistic was the conditional entropy h_{m+1}. Testing began with m = 0 and continued for increasing values of m until the null hypothesis (process is of order m) could not be rejected.

Results for two representative afferents are shown in Figure10, which depicts the h_{m} for afferent (*solid line*) and mean h_{m} for 49 realizations of the surrogate data (*dashed line*). For the unit shown in Figure10*A*, testing showed that afferent data had significantly smaller h_{m}(*) than surrogate data (□) for orders up to 3. For orders 4 and 5, there was no significant difference in the h_{m} (○). It was concluded that this afferent could be described by a third-order Markov process. Testing terminated for only 5 of 52 afferents with the null hypothesis being accepted. Testing also terminated if there was not sufficient data for the afferent. This is shown in Figure 10*B*, where the afferent and surrogate data sets were significantly different up to fifth-order. There was not sufficient data to test for higher orders, and thus m = 5 is only a lower-bound on the order of the process for this afferent. Testing terminated for the majority of afferents (47 of 52) in this way.

Table 1 summarizes the test for Markov order. The total number of afferents for which the order was at least m are listed for each m in the top row. Of these, the number of afferents for which the order was exactly m are shown in the bottom row. For example, all 52 afferents were at least order 2 (top row, m = 1, 2), with one afferent being exactly order 2 (bottom row, m = 2). Similarly, 51 were at least order 3, of which three afferents were exactly order 3. Thus, discarding the three afferents that were exactly order 3, only 48 afferents were tested for m = 4. Four of these could not be tested because of insufficient data, leaving 44 afferents that were at least order 4 (top row, m = 4). Of these afferents, only one was exactly order 4 (bottom row, m = 4). The test proceeded in this manner until afferents were exhausted. Most afferents (44 of 52) were fourth-order Markov or greater, and approximately a quarter of the units (12 of 52) were seventh-order Markov or greater.

### Count distributions

For nonrenewal processes, spike count distributions can provide additional information about the process and can complement higher-order interval analysis. Spike count distributions were analyzed from the count sequence defined by Equation 2. The SD ς_{C}, coefficient of variation CV_{C}, and variance-to-mean ratio F_{C} (also called Fano factor) of the spike count distributions were computed as a function of window duration T. T plays the same role in count analysis as interval order k in interval analysis.

Figure 11*A–C* shows the statistical measures for the same representative afferent fiber, whose interval analysis was described earlier (Fig. 3*A–C*). As with interval analysis, spike count distributions exhibited a minimum variance-to-mean ratio (Fano factor F_{C}) at an optimum count window of duration T_{min}. For the fiber shown in Figure 11*C*, T_{min} = 280 EOD periods (368 msec) with F_{C}(T_{min}) = 0.01 spikes. Mean spike count at T_{min} was 94 ± 0.98 spikes, i.e., a SD that is ∼1% of the mean (Fig.11*A,B*). Beyond this minimum T, counts became more irregular with increasing T because of the rapid increase in SD (Fig. 11*A*). The trends in ς_{C}, CV_{C}, and F_{C} were observed in all afferents and closely resemble those seen in interval analysis (Fig. 3). Figure11*D–F* summarizes the median value of ς_{C}, CV_{C}, and F_{C} for the population. The distribution of T_{min} across the population is shown in the histogram overlaying the Fano curves in Figure 11*F.*Mean T_{min} for the population was 233 ± 172 EOD periods (270 ± 200 msec), and mean F_{C}(T_{min}) for the population was 0.025 ± 0.022. At T_{min} when fibers demonstrate most regular counts, the distribution of CV_{C}(T_{min}) for the population is shown in Figure 4*B.* Twenty seven of 52 fibers had SDs that were <2% of mean spike count. Mean CV_{C}(T_{min}) for the population was 0.026 ± 0.023.

As was the case for interval analysis, the spike count analysis indicates that the afferents are significantly more regular than the three model spike trains (B, M_{0}, and M_{1}). The improvement in regularity measured as a ratio of the Fano factor of the model to that of the afferent at T_{min} is shown in Figure12. The figure is similar to Figure 6. Afferents had Fano factors that were on the average smaller than those of the models by 19.6 ± 10.5 (B, Fig.12*A*), 6.0 ± 4.7 (M_{0}, Fig. 12*B*), and 1.7 ± 1.0 (M_{1}, Fig.12*C*). Thus, the count analysis agrees with interval analysis in that the regularity of afferents cannot be explained by renewal models (B and M_{0}) or even a nonrenewal model that incorporates adjacent interval dependencies (M_{1}).

Although interval and count analysis are in agreement on the general trends seen in afferent data, it should be noted that there are differences between the two. The most significant difference lies in a determination of the time scale on which afferents may be considered most regular. From interval analysis, the mean time estimated for the population at their k_{min} was 176 msec (see Interspike interval analysis), which is smaller than the mean T_{min} estimated at 270 msec from the count analysis (see above).

### Weak signal detection

The regularity of primary electrosensory afferent spike trains is most pronounced for count windows on the order of ∼200 EOD periods (∼200 msec). Here, we illustrate the impact that this spike count regularity could have on detection performance, if the nervous system were to implement a detection algorithm using a comparable integration time constant. Full analysis of the optimal detection algorithm for the actual spatiotemporal profiles of naturally occurring electrosensory signals is beyond the scope of this paper. Here, we demonstrate the magnitude of the effect by considering a simplified signal detection experiment in which the stimulus is assumed to transiently increase the spike count on an individual afferent by a fixed number of spikes n_{s} within a fixed time window.

The detection algorithm is based on an ideal observer paradigm (Green and Swets, 1966) using spike count distributions (see Materials and Methods). It implements a binary hypothesis testing procedure illustrated schematically in Figure13*A.* In the absence of a stimulus, the spike count is drawn randomly from the baseline spike count distribution (Fig. 13*A*, *Baseline*). In the presence of the stimulus, the spike count is assumed to be drawn from a shifted version of the same distribution, with an offset of n_{s} spikes (Fig. 13*A*,*Baseline + Signal*). Given an observed spike count, the detection experiment decides which of the two distributions is most likely to have generated the observation. If the count exceeds a fixed threshold (Fig. 13*A*, *Threshold*), then we decide that a stimulus event has occurred. This could be a false alarm because of a baseline fluctuation being erroneously classified as a stimulus event (area marked as P_{fa} under the*Baseline* curve) or a correct detection if a stimulus was in fact present (area marked as P_{d} under the*Baseline + Signal* curve).

Simulated signal detection experiments were performed using afferent data from a representative unit, as well as three matched spike train models (B, M_{0}, and M_{1}). The results are shown in Figure13*B.* In all cases the detection probability P_{d} increased with increasing signal strength n_{s}, but it increased more rapidly with increasing number of added spikes for the afferent than for any of the three models. If we arbitrarily select P_{d} = 0.90 for comparison purposes (Fig. 13*B*,*horizontal dashed line*), we find that an addition of two or three spikes results in a 90% detection probability for the afferent spike train, whereas a comparable level of performance requires approximately five spikes for M_{1}, 10 spikes for M_{0}, and 18 spikes for the binomial process B. Specifically, comparing results for the afferent data and the best-matched renewal process model M_{0}, we see that the afferent spike train permits efficient and reliable detection of signals that are a factor of 3–5 times weaker than could be detected if baseline activity were generated by a renewal process.

## DISCUSSION

The principal finding of this study is the remarkable regularity of P-type afferents on intermediate time scales, given the high variability of their interspike intervals on short time scales. We measured variability of intervals using SD (ς_{I}), coefficient of variation (CV_{I}), and variance-to-mean ratio (F_{I}). Across the population of afferents the coefficient of variation of the first-order ISI, CV_{I}(1) averaged ∼0.44. Higher-order intervals I_{k}(which are the sums of k successive ISIs) showed a rapid decrease in CV_{I}. On intermediate time scales (k ≈ 50), CV_{I}(k) was <0.02 for most afferents (Fig. 3*E*). In contrast if the intervals were generated by a renewal process, CV_{I} would decrease as k^{−1/2}, and the expected CV_{I}on intermediate time scale would have been ∼3–5 times larger.

Although CV_{I} is a normalized measure of variability, it is easier to understand the large improvements in regularity of afferents by examining how the standard deviation of intervals ς_{I} changes with k (Fig.3*A*). For a renewal process ς_{I}increases as k^{1/2}, but for afferents, there is little increase in ς_{I} up to k ≈ 50 (Fig. 3*A*). It is as if the spike-generating process were keeping a check on the cumulative deviation of the successive ISIs from the mean ISI. Such a regularizing process could cancel deviations so that overall variability of k successive intervals is approximately the same as the variability of a single interval. If such a process were responsible for the observed regularity it would appear that the process is able to maintain this regularity only over intermediate time scales. This is seen most clearly in the variance-to-mean ratio curve F_{I}(Fig. 3*C*). The F_{I} curve decreases sharply with increasing k, whereas a renewal process has constant F_{I}. For the afferent data, F_{I} decreases until it reaches a minimum (for some afferent specific k_{min}), after which it loses regularity and increases with increasing k. The time scale on which the spike regularity is most pronounced relative to a renewal process corresponds to interval order k_{min}. The distribution of k_{min} in the population of afferents is shown in the *gray histogram* of Figure3*F.* The population mean was
_{min} = 42, which corresponds to ∼175 msec. Spike count distributions (Fig. 11) demonstrated similar trends, although the estimated time scale over which spike counts were most regular was somewhat larger (270 msec).

The rapid decrease in variability has implications for signal detection. For a renewal process, to achieve a similar reduction in variability would require a time scale of about k_{min}^{2}, which is ∼1700 ISIs, or several seconds. Therefore, in situations where sensory signals are extremely weak, the regularization process allows reliable detection using smaller integration times.

### Interval and count variability in other sensory systems

The coefficient of variation of the ISI distribution has long been the standard measure for specifying spike train variability. For the benchmark Poisson process, CV_{I} = 1 for all interval orders. Spike trains with CV_{I} > 1 (< 1) are more irregular (regular) than a Poisson process. Spontaneously active neurons demonstrate a wide range in ISI variability with reported values ranging over two orders of magnitude from 0.02 (Ratliff et al., 1968) to ∼3 (Teich et al., 1997). In comparison, the values we report here for P-type electroreceptor afferents had a mean CV_{I}(1) of 0.44 and ranged from 0.15 to 0.79, reflecting moderate to high short-term variability.

For renewal processes, the ISI distribution is sufficient to completely describe variability on all time scales (Cox, 1962). For nonrenewal spike trains, such as P-type afferents reported here, it is necessary to examine higher-order distributions. Towards this end, interval distributions I_{k} or count distributions C_{T} and their statistical properties can provide useful information about variability. Higher-order interspike interval distributions were first introduced by Rodieck et al. (1962), but with the exception of second-order interval distributions I_{2} (Teich and Khanna, 1985) they have not been widely used. This is unfortunate because higher-order interval distributions convey useful information about variability on multiple time scales. On the other hand, spike count distributions have been more widely used in recent years. However, they have been used primarily to characterize spike trains for which the Fano factor exhibits power-law growth with count window T, on time scales extending beyond tens of seconds (for review, see Teich et al., 1996). Such spike trains exhibit a high degree of irregularity (F_{C} > 1) that increases with counting time. On the time scales that are of interest to us (0.1–1 sec), most of the reported work on spike count distributions are restricted to driven responses (Teich and Khanna, 1985; Young and Barta, 1986;Relkin and Pelli, 1987; Shofner and Dye, 1989; Softky and Koch, 1993; Baddeley et al., 1997; Shadlen and Newsome, 1998) and cannot be directly compared with the results for baseline responses.

### History effects and nonrenewal nature of afferent spike trains

Most experimentally observed spike trains exhibit correlations and memory effects in the ISI sequence. That is, they are described by nonrenewal processes. This is true for all P-type afferents reported here. All afferents had adjacent ISIs that were strongly anticorrelated. Anticorrelations have been noted in spike trains of P-type units in other species of electric fish such as the gymnotid*Sternopygus* (Bullock and Chichibu, 1965) and*Steatogenes* (Hagiwara and Morita, 1963). They have also been reported in many other sensory systems including visual (Kuffler et al., 1957), auditory (Johnson et al., 1986; Lowen and Teich, 1992), and somatosensory (Amassian et al., 1964), and in motor systems (Calvin and Stevens, 1968). Whereas anticorrelation improves regularity, our results suggest that it is not sufficient to account for the dramatic improvements in regularity observed in P-type afferents. We demonstrated that a first-order Markov process M_{1} that incorporated the anticorrelations observed in the afferent data could not reproduce the reduction in variability demonstrated by the afferent spike trains with increasing T (Fig. 11). This suggested that dependencies between ISIs persisted over many intervals.

To assess the duration of history-dependent effects, serial correlation coefficients were evaluated. However, SCCs suffer from two problems. First, even a first-order Markov process can have statistically significant correlations that persist over many intervals (Fig. 8). Second, even if a process is higher-order Markov (Fig. 9*B*), the use of serial correlation coefficients can average out effects of dependencies for lags smaller than the Markov order (Fig.9*A*). Thus, Figures 8 and 9 demonstrate that the SCCs cannot be used to determine the order of the Markov chain, and an explicit test of Markov order is required (Nakahama et al., 1972). When explicitly tested for the order of the underlying Markov chain, the majority of afferents were typically fourth-order or greater (as also reported by van der Heyden et al., 1998).

The high degree of regularity exhibited by afferents is a consequence of a nonrenewal process that keeps a check on the deviation from the mean firing rate over many ISIs. Typically, when an interval longer than the mean ISI occurs (a “credit”), the next interval will be shorter than the mean (a “debit”). This gives rise to the strong long–short anticorrelations observed in the data. However, when the credit does not get paid back immediately, it is not forgotten. Rather, it is paid back eventually. This is seen from Figure 9*B*where we evaluated the distribution of Y given the previous occurrence of X in the ISI sequence (X, 2, 2, 2, Y). The joint distribution shows that the long–short distribution of intervals persists even when there are a number of intervening intervals where the debits and credits do not cancel each other. It is likely that the demands placed on a physiological system for maintaining precision or regularity of spiking may be more easily achieved by making adjustments in timing over many intervals rather than over a few intervals. At present the mechanism causing this remarkable degree of precision is unknown. Although mathematical models of the spike-generating process in P-type afferents have been constructed (Longtin and Racicot, 1997b; Longtin, 1998), they do not explain the long-range interval correlations noted here.

### Implications for detection of weak sensory signals

The increased spike train regularity that we observe for P-type afferents on intermediate time scales has important implications for the ability of the fish to detect weak electrolocation targets. P-type afferents encode amplitude modulations of the local transdermal voltage caused by nearby objects (Scheich et al., 1973;Bastian, 1981) (for review, see Zakon, 1986). Objects in the water near the fish modulate the transdermal potential and thus modulate the per-cycle firing probability of P-type afferents. Objects that cause a large change in firing probability are easily detected by the animal. However, when the object is small, or far away, or has low electrical contrast with the surrounding water, then the resulting small change in per-cycle firing probability can be obscured by statistical fluctuations in the baseline spike activity. Early studies have established that the behavioral threshold for *Apteronotus* is <1 μV/cm (Knudsen, 1974). Bastian (1981) extrapolated his data to show that in this range, the change in firing rate of P-type afferents is ∼1 spike/sec. Thus, at threshold levels of performance the expected change in firing rate is <1% (based on a mean discharge rate of 294 spikes/sec as reported by Bastian, 1981).

In behavioral studies of prey capture performed in our laboratory, we have shown that *Apteronotus* can detect small water fleas (*Daphnia magna*, 2–3 mm in length) at a distance of ∼2 cm from the electroreceptor array (Nelson and MacIver, 1999). At this distance, the spatial profile of local transdermal voltage change is a Gaussian-like bump with a full width at half maximum of ∼2 cm (Rasnow, 1996; Nelson and MacIver, 1999). During prey capture behavior the fish is typically moving with a relative velocity of ∼10 cm/sec relative to the prey. Thus, as the electrosensory image passes over the receptor array, an individual electroreceptor organ experiences a transient change in transdermal voltage that lasts ∼200 msec. Interestingly, this temporal duration is well matched to the time scale on which P-type afferents show the greatest increase in spike train regularity.

Based on computer reconstructions of electrosensory images from our behavioral studies, we estimate that the peak change in per-cycle firing probability of P-type afferents at the time of prey detection is at most a few percentage of the baseline probability (Nelson and MacIver, 1999). Assuming a typical baseline spike rate of 300 spikes/sec, this corresponds to approximately one extra spike of a total of 60 spikes expected in a 200 msec time window on an individual afferent. This is a weak signal and one that could potentially be difficult to detect given the relatively high variability observed in the first-order ISI distribution.

The neural computations required to detect the dynamically changing spatiotemporal profile across a population of electrosensory afferents could plausibly be implemented by circuitry in the brainstem electrosensory nucleus, the electrosensory lateral line lobe (Shumway, 1989; Metzner et al., 1998). A full analysis of the optimal detection algorithm is beyond the scope of this paper. However, in performing such an analysis, it will be extremely important to take into account the intermediate-term regularity of P-type afferent spike trains afforded by the underlying nonrenewal process. If we performed such an analysis assuming that P-type afferents could be adequately modeled by matching the first-order ISI distribution and assuming a renewal process model, we would overestimate the behavioral threshold by a large factor (Fig.13*B*). Thus, it is important for sensory physiologists and neural modelers to consider potential effects of nonrenewal statistics when analyzing detection thresholds for weak sensory signals in noisy spike trains.

## Appendix

### Spike count distributions for discrete-time renewal processes

For spike trains arising from a discrete-time renewal process, the distribution of spike counts can be calculated from the distribution of interspike intervals. To illustrate this, we first present results for the binomial spike train and then derive expressions relating count and interval distributions for arbitrary renewal processes.

For a binomial process with constant per-cycle probability of firing p, the probability of observing an interval of j cycles between successive spikes is:
Equation 10which is the geometric distribution with mean
_{1} = 1/p. In general, the probability of observing an interval j between k successive spikes (i.e., the k th order interval distribution) is:
Equation 11with mean
_{k} = k/p. Equations 10and 11 are the discrete time analogs of the exponential and k th-order gamma distributions, respectively.

For a binomial process, the number of spikes in a count window of duration T cycles is distributed according to:
Equation 12where Ĉ_{T}(k) is the probability of finding k spikes in T cycles.

For an arbitrary discrete-time spike-generating process, if counting is performed so that the window starts on the cycle immediately after a spike, then the distribution for the number of spikes in a window of size T, N̂_{T}, is related to the distribution of the k th order spike intervals, S_{k}, by:
Equation 13In other words, the probability that there are fewer than k spikes in a count window of size T is the same as the probability that the k th order interval is greater than T. The count distribution and interval distribution are therefore related by:
Equation 14For a renewal process, the ISIs are independent and identically distributed with some known probability density function I_{1}(j). Because any k th order interval of length r can be expressed as r = ∑_{i=1}^{k} j_{i}, where the j_{i} are the intervening ISIs, the renewality condition implies that the k th order interval distribution is the k -fold convolution of I_{1} with itself (Feller, 1957). That is,
Equation 15For renewal spike trains, we can use Equation 15 to express I_{k+1} in terms of I_{k}. Introducing this into Equation 14, and after some routine algebra, we have:
Equation 16From Equation 15 it can be seen that I_{k} is completely specified if I_{1} is known, and hence, from Equation 16 count distributions are also known. That is, for a renewal process the ISI distribution I_{1} is sufficient to completely describe the count distributions.

As an example, we can obtain the distribution of counts for a binomial spike train (Eq. 12) from the ISI distribution (Eq. 10). Repeated application of Equation 15 with I_{1} the geometric density given by Equation 10 yields the k th order interval distribution given by Equation 11. Furthermore, inserting Equations 10and 11 into Equation 16 yields the binomial distribution given by Equation 12.

For a nonrenewal process, the chief difficulty lies in relating the ISI distribution I_{1} to the k th-order interval distributions. That is, neither the convolution expression given by Equation 15 nor the expression for count distribution given by Equation 16 are applicable. Therefore, Equation 14 can be evaluated only if I_{k} for all k ≥ 1, are known. Thus, for a nonrenewal process, the ISI distribution does not provide information about counts or intervals on multiple time scales.

## Appendix

### Measures of variability for renewal processes

Let I_{1} be the ISI distribution for a renewal process, and let
_{1} and Var (I_{1}) be its mean and variance, respectively. Then, as the k th-order interval is the sum of k independent random variables, the mean and variance of I_{k} are k
_{1} and k Var (I_{1}), respectively. Thus, the coefficient of variation CV_{I}(k) is:
Equation 17That is, the coefficient of variation decreases as k^{−1/2}. Likewise, the variance-to-mean ratio F_{I}(k) is:
Equation 18For the binomial model, it follows from Equation 11 that
_{k} = k/p and Var (I_{k}) = k(1 − p)/p^{2}. Therefore, CV_{I}(k) =
and F_{I}(k) = (1 − p)/p.

For a nonrenewal process, Var (I_{k}) can be related to Var (I_{1}) and the correlation coefficients ρ_{l} by Var (I_{k}) = k Var (I_{1}){1 + 2∑_{l=1}^{k−1}(1 − (
))ρ_{l}}.

The coefficient of variation and Fano factor expressions for spike count distributions are more complicated than those for the intervals, with the exception of binomial (and Poisson) processes for which they are easily calculated. For the binomial spike train, it follows from Equation 12 that
_{T} = Tp and Var (C_{T}) = Tp (1 − p). Thus, CV_{C}(T) =
, and falls as T^{−1/2}. The Fano factor F_{C}(T) = 1 − p, and is constant for all T. Any process with F_{C}(T) smaller than this value will exhibit greater regularity in spike counts than the binomial process, and conversely, more irregular processes will have larger values.

For small T, F_{C}(T) for any discrete time process will tend to 1 − p where p is the probability of firing in the sampling interval. That is, for counting times T → 1, the process is nearly binomial. For a continuous time process, when the count window Δt → 0 then p → 0. Thus, F_{C}(T) → 1, which is the Fano factor for a Poisson process. That is, we expect spike counts in small count windows to be as irregular as a Poisson process. For renewal processes, the Fano factor for large T asymptotically approaches a constant value that is related to the coefficient of variation of the ISI by the relation, F_{C}(T) ≈ CV_{I}(1)^{2}(Cox and Lewis, 1966). That is, the variability of spike counts is attributable to the variability in interspike intervals.

## Footnotes

This research was supported by National Institute of Mental Health Grant R01 MH-49242. We thank Noura Sharabash and Zhian Xu for experimental support and Josina Goense for helpful comments and suggestions on this manuscript.

Correspondence should be addressed to Rama Ratnam, 3317 Beckman Institute, 405 North Mathews Avenue, Urbana, IL 61801. E-mail:ratnam{at}uiuc.edu.