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 (Ptype) electrosensory afferents encode amplitude modulations of the fish's selfgenerated electric field and provide information necessary for electrolocation. This study characterizes the statistical properties of baseline spike activity in Ptype afferents of the brown ghost knifefish, Apteronotus leptorhynchus. Shortterm 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 highorder interval analysis, count analysis, and Markovorder 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 fourthorder 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 firstorder 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 firstorder ISI distribution only provides information about shortterm variability on time scales comparable to the mean ISI. Longterm variability, over time periods containing multiple spikes, must be measured using other techniques, such as analysis of higherorder 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 longterm 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, higherorder interval and count distributions can be derived knowing only the firstorder 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 historydependent effects in the ISI sequence. In such cases, the firstorder ISI distribution does not provide sufficient information to predict longterm spike train variability or to set limits on signal detection performance.
In this paper, we analyze the variability of baseline spike activity recorded from Ptype (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 selfgenerated 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). Ptype 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 inApteronotus 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 Ptype afferents by only ∼1% (Nelson and MacIver, 1999). In this study we show that Ptype afferent spike trains are irregular as judged by the firstorder 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 Ptype 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 (MS222; Sigma, St. Louis, MO) and immobilized with an intramuscular injection of 3 μl 10% gallamine triethiodide (Flaxedil; Sigma). Ptype 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 timestamped 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 quasisinusoidal EOD waveform with a fundamental frequency f that ranges from 750 to 1000 Hz depending on the individual. Ptype units fire at most once per EOD cycle and randomly skip cycles between successive spikes. On average, a typical unit fires on about onethird of the EOD cycles. This ratio is referred to as the percycle probability of firing p. In the presence of a stimulus, the percycle 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 discretetime 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 percycle 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)
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 andD, 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:
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 Ptype afferent is that it is a “probability coder.” Namely, the afferent fires irregularly with a percycle probability p that is constant under baseline conditions, but is subject to modulation in the presence of external stimuli. Thus, the simplest model of Punit baseline activity is a process that emits a spike in any given EOD cycle with constant probability p independent of previous history. In a continuoustime 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 discretetime 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 Poissondistributed 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. Percycle firing probability p remains unchanged because shuffling preserves the total number of spikes and EOD periods.
Zerothorder 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 zerothorder Markov process.
Firstorder 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 firstorder Markov process and will be denoted by M_{1}.
The binomial (B) and zerothorder 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 firstorder 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:
Another useful measure is the variancetomean ratio of I_{k}, denoted F_{I}(k):
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 variancetomean 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:
For spike count distributions, the variancetomean ratio is called the Fano factor (Fano, 1947). It is denoted F_{C}(T) and defined as:
Correlation analysis
Historydependent effects were analyzed by considering the serial correlation coefficient (SCC) ρ_{l} of the firstorder sequence S_{1}, where l is the lag in terms of the number of intervening intervals. The ρ_{l} were computed from:
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. Higherorder historydependent 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 zeroorder 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 firstorder 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:
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 onesided, thep value for the test is p = r/(R_{s} + 1). The null hypothesis was rejected ifp ≤ 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 lowerbound 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 Ptype 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 Ptype 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 (
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 firstorder 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 higherorder 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 variancetomean 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
The trends observed in ς_{I}, CV_{I} and F_{I} for this representative unit were consistent across the entire population. Figure 3
D–Fillustrates 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 variancetomean 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 (
The higherorder 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 (percycle 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 andA2 (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
Nevertheless, the binomial model serves as a useful benchmark because the coefficient of variation CV_{I} and variancetomean ratio F_{I} can be computed analytically (see Appendix ). For the k thorder interval distribution, ς_{I}(k) =
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 variancetomean ratios that were on the average 50 times smaller (mean, 50 ± 27) than the F_{I} for binomial spike trains.
Zerothorder 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 zerothorder 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 variancetomean ratio F_{I}can be computed analytically (see Appendix ): CV_{I}(k) = CV_{I}(1)/
Although M_{0} is more regular than the binomial model B, it still does not capture the intermediateterm regularity of the afferent data. At the minimum of the variancetomean 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 Ptype 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 andC2.
The simplest nonrenewal model that incorporates the observed long–short dependency is the firstorder 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 secondorder 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 (
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 firstorder Markov model M_{1} (dashdot 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 firstorder (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 zerothorder Markov process for the same afferent as shown in Figure5 C2. The zerothorder 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, longterm 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 thirdorder 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 fifthorder. There was not sufficient data to test for higher orders, and thus m = 5 is only a lowerbound 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 fourthorder Markov or greater, and approximately a quarter of the units (12 of 52) were seventhorder Markov or greater.
Count distributions
For nonrenewal processes, spike count distributions can provide additional information about the process and can complement higherorder interval analysis. Spike count distributions were analyzed from the count sequence defined by Equation 2. The SD ς_{C}, coefficient of variation CV_{C}, and variancetomean 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 variancetomean 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 theBaseline curve) or a correct detection if a stimulus was in fact present (area marked as P_{d} under theBaseline + 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 bestmatched 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 Ptype 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 variancetomean ratio (F_{I}). Across the population of afferents the coefficient of variation of the firstorder ISI, CV_{I}(1) averaged ∼0.44. Higherorder 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 spikegenerating 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 variancetomean 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
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 Ptype electroreceptor afferents had a mean CV_{I}(1) of 0.44 and ranged from 0.15 to 0.79, reflecting moderate to high shortterm 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 Ptype afferents reported here, it is necessary to examine higherorder distributions. Towards this end, interval distributions I_{k} or count distributions C_{T} and their statistical properties can provide useful information about variability. Higherorder interspike interval distributions were first introduced by Rodieck et al. (1962), but with the exception of secondorder interval distributions I_{2} (Teich and Khanna, 1985) they have not been widely used. This is unfortunate because higherorder 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 powerlaw 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 Ptype afferents reported here. All afferents had adjacent ISIs that were strongly anticorrelated. Anticorrelations have been noted in spike trains of Ptype units in other species of electric fish such as the gymnotidSternopygus (Bullock and Chichibu, 1965) andSteatogenes (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 Ptype afferents. We demonstrated that a firstorder 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 historydependent effects, serial correlation coefficients were evaluated. However, SCCs suffer from two problems. First, even a firstorder Markov process can have statistically significant correlations that persist over many intervals (Fig. 8). Second, even if a process is higherorder 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 fourthorder 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 Bwhere 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 spikegenerating process in Ptype afferents have been constructed (Longtin and Racicot, 1997b; Longtin, 1998), they do not explain the longrange interval correlations noted here.
Implications for detection of weak sensory signals
The increased spike train regularity that we observe for Ptype afferents on intermediate time scales has important implications for the ability of the fish to detect weak electrolocation targets. Ptype 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 percycle firing probability of Ptype 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 percycle 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 Ptype 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 Gaussianlike 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 Ptype 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 percycle firing probability of Ptype 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 firstorder 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 intermediateterm regularity of Ptype afferent spike trains afforded by the underlying nonrenewal process. If we performed such an analysis assuming that Ptype afferents could be adequately modeled by matching the firstorder 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 discretetime renewal processes
For spike trains arising from a discretetime 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 percycle probability of firing p, the probability of observing an interval of j cycles between successive spikes is:
For a binomial process, the number of spikes in a count window of duration T cycles is distributed according to:
For an arbitrary discretetime spikegenerating 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:
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 thorder 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
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 − (
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
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 MH49242. 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. Email:ratnam{at}uiuc.edu.