## Abstract

In vertebrate auditory systems, the conversion from graded receptor potentials across the hair-cell membrane into stochastic spike trains of the auditory nerve (AN) fibers is performed by ribbon synapses. The statistics underlying this process constrain auditory coding but are not precisely known. Here, we examine the distributions of interspike intervals (ISIs) from spontaneous activity of AN fibers of the barn owl (*Tyto alba*), a nocturnal avian predator whose auditory system is specialized for precise temporal coding. The spontaneous activity of AN fibers, with the exception of those showing preferred intervals, is commonly thought to result from excitatory events generated by a homogeneous Poisson point process, which lead to spikes unless the fiber is refractory. We show that the ISI distributions in the owl are better explained as resulting from the action of a brief refractory period (∼0.5 ms) on excitatory events generated by a homogeneous stochastic process where the distribution of interevent intervals is a mixture of an exponential and a gamma distribution with shape factor 2, both with the same scaling parameter. The same model was previously shown to apply to AN fibers in the cat. However, the mean proportions of exponentially versus gamma-distributed intervals in the mixture were different for cat and owl. Furthermore, those proportions were constant across fibers in the cat, whereas they covaried with mean spontaneous rate and with characteristic frequency in the owl. We hypothesize that in birds, unlike in mammals, more than one ribbon may provide excitation to most fibers, accounting for the different proportions, and that variation in the number of ribbons may underlie the variation in the proportions.

## INTRODUCTION

The central auditory system, like other sensory systems, encodes stimulus information in the timing and/or rate of action potentials (spikes). In vertebrates, the translation of a graded potential across the hair cell receptor membrane into a stochastic train of spikes in the primary afferents, the auditory nerve (AN) fibers, is performed by the first synapse in the auditory system, a ribbon synapse specialized for tonic release of neurotransmitter at high rates (for reviews, see Fuchs 2005; Moser et al. 2006a, b; Nouvian et al. 2006; Prescott and Zenisek 2005; Sterling and Matthews 2005). A better understanding of this process will promote a more thorough understanding of the way in which AN fibers encode acoustic stimuli and of the constraints of temporal coding. One approach by which useful information about this translational process can be derived is to analyze the detailed timing of AN fiber spikes and to develop models that account for the observations. Although the approach is rather indirect, because only the final output of a cascade of intervening steps is scrutinized, the hair cells and synapses can operate in an undisturbed environment while the measurements are made.

The spiking activity of AN fibers is commonly modeled as a stochastic point process (Perkel et al. 1967). Frequently, an excitatory stochastic process is assumed to generate excitatory events that trigger spikes unless the AN fiber is refractory. In the case of spontaneous activity, where in some AN fibers of mammals and birds discharge rates can exceed 100 spikes/s (Gummer 1991b; Köppl 1997a; Liberman and Kiang 1978; Richter et al. 1995; Taberner and Liberman 2005; Winter and Palmer 1991), and with the notable exception of avian AN fibers showing preferred intervals, the excitation is generally assumed to be provided by a homogeneous Poisson point process (for reviews, see Delgutte 1996; Mountain and Hubbard 1996). To account for rate trends and phenomena seen over longer time scales (>100 ms out to tens of seconds), others have assumed a fractal doubly stochastic Poisson point process (Jackson and Carney 2005; Lowen and Teich 1991, 1992). Both types of process lead to identical interspike interval (ISI) distributions (Teich et al. 1990). The observed deviations of measured ISI distributions from exponential distributions are therefore often attributed to the refractory properties of AN fibers (Carney 1993; Gaumond et al. 1983; Kiang et al. 1965; Li and Young 1993; Miller and Wang 1993; Young and Barta 1986; Zhang et al. 2001). However, we have recently shown that the ISI distributions obtained from spontaneous activity of cat AN fibers are not compatible with these assumptions (Heil et al. 2007). We also showed that, instead, a simple and physiologically plausible modification of an excitatory Poisson process, in combination with refractoriness that is consistent with experimental data (Brown 1994; Cartee et al. 2000; Dynes 1996; Miller et al. 2001; Morsnowski et al. 2006; Shepherd et al. 2004), provided a description that was nearly 2 orders of magnitude better than that achieved by a Poisson process. Furthermore, the description could fully account for the observed ISI distributions. The modification proposed by us still suggests a homogeneous excitatory stochastic process, but one where the distribution of intervals between excitatory events is a mixture of an exponential and a gamma distribution with shape factor 2, both with the same scaling parameter (or rate). This mixture can emerge from a single primary homogeneous Poisson point process, if some of its primary events are subsequently eliminated, obeying a “never-two-in-a-row” rule, and do not contribute to the ultimate excitatory events (Heil et al. 2007). As likely physiological correlates of such excitatory events, we suggested transmitter release events from the hair cell that give rise to excitatory postsynaptic potentials of amplitudes sufficient to exceed the threshold for spike initiation in the afferent fiber unless the fiber is refractory.

Here, we ask whether our findings of a mixed distribution of intervals between excitatory events in cat AN fibers might apply more broadly. We focus on birds, another group of vertebrates with highly developed auditory systems. Tall hair cells in the basilar papilla of birds are thought to functionally correspond to inner hair cells of mammals (review in Manley and Köppl 1998). Their synaptic ribbons, also called presynaptic dense bodies or synaptic balls in early papers, are very similar in appearance to those of mammalian inner hair cells (Chandler 1984; Fischer 1994; Martinez-Dunst et al. 1997; Takasaka and Smith 1971; Tanaka and Smith 1978). We chose to study the discharge behavior of AN fibers in the barn owl (*Tyto alba*). The barn owl is a well-studied species with regard to its auditory system, because of its nocturnal hunting habits and the associated specializations for sound localization. At the cochlear level, the barn owl shows an extended representation of high frequencies (4–10 kHz), associated with an unusually high number of ∼31,000 afferent nerve fibers and extreme temporal precision in the form of high-frequency phase locking (Köppl 1997b,c; Köppl et al. 1993). Finally, long records of spontaneous activity of well-characterized auditory primary afferents were already available (Köppl 1997a). Taking together, this made the barn owl the bird of choice for the purpose of this study.

## METHODS

The data for this study are a subset of those published by Köppl (1997a), where details of the surgical procedures, of the acoustic stimulation, and of recording techniques can be found. In brief, seven adult barn owls (*Tyto alba guttata*, 4 females and 3 males) were used. General anesthesia was induced by intramuscular injections of 3 mg/kg xylazine (Rompun) and 4 mg/kg ketamine hydrochloride and maintained by additional doses as required; the ECG and breathing cycle were continuously monitored via recordings from needle electrodes in wing and leg muscles. At 6- to 12-h intervals, a dose of diazepam (valium, 0.8 mg/kg) was given instead of xylazine. Cloacal temperature was kept at 39–40°C with the aid of a feedback-controlled heating pad wrapped around the owl's body. A metal pin was glued to the skull for head fixation, and the right auditory brain stem was exposed by opening the skull and meninges and aspirating part of the cerebellum.

### Stimulation and recording procedures

For data acquisition, the owls were situated in a custom-built double-walled sound-attenuating chamber. Closed sound systems were sealed to both ears and individually calibrated for each animal. Single AN fibers were recorded with glass micropipettes filled with 3 M KCl or 2 M NaCl, with resistances mostly between 50 and 100 MΩ. Electrodes were positioned just above the surface of the brain stem under visual control and remotely advanced via a hydraulic microdrive and a piezo-stepper. Recording signals from a WPI 767 electrometer could be variably filtered, and action potentials were fed via a threshold discriminator (Schmitt-trigger) to a custom-built computer interface. Spike times were taken as the instances at which the amplified and filtered electrode signal crossed the Schmitt-trigger level and were stored on disc with 5-μs or, in the first experiment, with 10-μs resolution for off-line analysis.

The single-unit nature of the recordings was verified by the presence of a minimum interval between successive spikes and a uniform spike amplitude at any given time. With the surgical approach used, most recordings were from AN fibers, but output fibers from the nucleus magnocellularis and, in rare instances, of the second-order nucleus laminaris could also be encountered. A reliable classification scheme, based on a number of response properties, was developed that allowed the post hoc separation of these neuron classes (see Köppl 1997a). Only recordings classified as AN fibers were used in this study. In addition, the sample was further restricted to a subset of AN fibers for which extensive recordings of spontaneous activity, ≤6 min, were obtained.

### Data analysis

For off-line analysis, spike times were imported into Excel 2000 spreadsheets and Visual Basic for Applications (VBA) routines (Microsoft, Redmond, WA) were developed to compute the more complex measures and to simulate data, as explained below. The Solver module of Excel was used to fit the models to the data using the Newton procedure.

##### GENERATION OF SUBSAMPLES AND SELECTION OF PERIODS OF STABLE SPONTANEOUS RATE FROM WITHIN THE SAMPLES.

The spontaneous activity of AN fibers can be subject to some rate changes over longer time scales, reflected for example in large Fano factors (≫1) for long counting times (Lowen and Teich 1992; see also Jackson and Carney 2005). To minimize any potential bias of such rate changes on models and parameter estimates, we applied strict criteria to select from each sample one or more periods of sufficiently stable spontaneous rate, very similar to the criteria and the procedures developed in our previous study of cat AN fiber spontaneous activity (Heil et al. 2007). Because of the considerably longer periods of time over which the spontaneous activity of barn owl AN fibers was recorded (≤6 min; see above), application of these criteria would have led to extensive pruning of the samples and thus to a loss of large numbers of ISIs for analysis. For probing the models (described in results), we therefore decided to cut long samples, with >4,000 ISIs (*n* = 82), into shorter subsamples, comprising 4,000 ISIs each, and a final subsample comprising the remaining ISIs from the full sample. Shorter samples, with <4,000 ISIs (*n* = 55), were left intact, but for convenience will also be referred to as subsamples. Each subsample was subjected to the stability analysis, as described next. We first plotted for each subsample the cumulative sum of the ISIs as a function of ISI number, in the sequence in which they were recorded. Longer-term fluctuations in spike rate become obvious in such plots as smooth changes in the slope of the function, and the offsets of straight line fits can deviate substantially from zero. We fitted straight lines to such functions for each subsample and computed the RMS of the residuals for each such fit. We also separated each subsample of spontaneous activity into two portions, one covering the first and the other the second half of the ISIs of the complete subsample, and computed the mean ISI and the associated SD for each portion. A subsample of spontaneous activity had to fulfil two criteria to be considered stable. First, the difference between the means of the ISIs during the first and the second half of the subsample had to be smaller than twice their probable difference, *D*, given by *D* = (SD_{1}^{2}/*n*_{1}+SD_{2}^{2}/*n*_{2})^{0.5}, where SD_{1} and SD_{2} are the SD from the first and second half of the subsample, respectively, and *n*_{1} and *n*_{2} are the corresponding numbers of ISIs. The difference between *n*_{1} and *n*_{2} was either one or zero. For *n*_{1} and *n*_{2} greater than ∼30, the means can be expected to be normally distributed, so this criterion has an error probability of ∼5%. In other words, ∼5% of random spike trains with stable longer-term rate can be expected to not meet this criterion. Second, the value of the offset (in ms) of the straight-line fit to the function relating the cumulative ISI to ISI number had to be smaller than 30 times the value of the square root of the RMS (also in ms) computed from the fit to all ISIs. If one or both of the two criteria were not fulfilled by a given subsample, the subsample was pruned until the remaining period fulfilled both criteria. Pruning started by discarding the ISIs from the last to that at which the cumulative function last intersected the straight-line fit. If after this step, the criteria were still not fulfilled, pruning was continued but this time from the beginning of the sub-sample, discarding the ISIs from the first to that at which the cumulative function first intersected the straight-line fit. If necessary, pruning continued in this alternating fashion, until the remaining period fulfilled both criteria. If no such period existed, the entire subsample was discarded.

Only periods of subsamples that passed the stability criteria and that comprised a minimum of 400 ISIs were used to fit the models. Of the 402 subsamples, generated from the recordings of spontaneous activity of 137 AN fibers (see results), 262 (65.2%) were pruned in this way, and 399, from 136 fibers, retained >400 ISIs. On average, 16.7% of the ISIs were removed from the 399 subsamples by the pruning. Application of the same stability criteria to a set of simulated spike trains generated with a constant rate resulted in pruning of 65 of 399 (16.3%) such trains and an average removal of ∼2.4% of the ISIs.

From each of the 399 remaining subsamples, we obtained estimates of the model parameters. However, although the individual subsamples constitute independent observations, the parameter estimates obtained from the different subsamples of a given AN fiber are not independent. Therefore for estimating parameters across the AN fiber population, we first computed the 82 means, and the SD, of the parameter estimates derived from the different subsamples of a given AN fiber and combined these means with the estimates from the 54 uncut shorter samples, which had passed the stability criteria, to compute the average across all fibers.

##### MODEL FITTING TO THE CUMULATIVE DISTRIBUTION FUNCTION.

Cumulative distribution function (CDF), probability density function (PDF), and hazard function (HF) are mathematically equivalent (Perkel et al. 1967). With finite numbers of realizations, however, empirical PDFs or HFs require binning of ISIs and thus contain less information, and they are inherently noisier than the corresponding CDF (see Fig. 2). Therefore the models were fitted directly to the empirical CDF, henceforth termed sample CDF. Note that a model that fits a sample CDF must also fit the corresponding PDF and HF and that a model that does not fit a sample CDF also does not fit the corresponding PDF or HF, even though the mismatch may be difficult to see with the latter two functions because they are so noisy when based on typical numbers of recorded ISIs. A sample CDF was obtained by sorting all ISIs from a given subsample, or from the period of a subsample that had passed the rate stability criteria, by their magnitudes. The cumulative probability, *P*_{i}, of a given ISI was calculated using the margin correction *P*_{i} = *i*/(*n* + 1), where *i* is the rank of the *i*th ISI after sorting them by magnitude, and *n* is the total number of ISIs in the sample (Sachs 2002).

A given model was fitted to a sample CDF by means of the adaptive Newton's method, using a four-stage procedure. During the first three stages, we minimized the same cost function as in Heil et al. (2007). For each data point, we calculated the vertical difference (in units of probability) as well as the horizontal difference (in units of time) between the sample CDF and the model CDF. The cost function was the product of the sums of the squared vertical differences and the squared horizontal differences, the latter being weighted with the survival function. During the fourth and final stage, only the sum of the squared vertical differences was minimized, because it is the well-known, classical, measure of the goodness-of-fit (Stephens 1986).

We used a similar strategy as in Heil et al. (2007) to guide the fitting routine in its convergence on the absolute minimum of the cost function with a reasonable combination of parameter estimates. As explained in results, four free parameters are required to fit the full version of the gamma-mixture model. These parameters are an absolute refractory period, *t*_{D}, a time constant of a relative refractory period, 1/λ_{R}, a scaling parameter, or rate of the presumed excitatory process, λ_{E}, and a parameter, *b*, which captures the relative proportions of different components in a presumed mixture of distributions of excitatory events. During the first stage of the four-stage procedure, only *t*_{D} and 1/λ_{R} were free parameters, whereas *b* was held fixed and λ_{E} was estimated from the mean ISI by the method-of-moments (Ayyub 1997). The estimate was updated with every iterative step of the adaptive Newton procedure using this method. During the second stage, *t*_{D}, 1/λ_{R}, and λ_{E} were free parameters, whereas *b* was still fixed, and during the third and fourth stages, all parameters were free. For fits of the Poisson model, for which *b* = 0, the third stage was redundant. As explained in the results, we also fitted a reduced version of each model to the data. In these versions, the refractory period is approximated by a single parameter, *t*_{RP}. In the fitting procedure, this parameter replaces *t*_{D}, whereas 1/λ_{R} was set to, and kept at, zero.

The initial starting values for parameters *t*_{D} or *t*_{RP} were set to 90% of the shortest ISI in the sample and the starting values for 1/λ_{R} and for *b* were set to 0 so that the starting conditions for the gamma-mixture model were identical to those for the Poisson model. The starting value for the scaling parameter was calculated by the methods-of-moments from these values and from the sample's mean ISI (see *Eqs. 5* and *9*). Extensive fits of simulated random data generated according to the gamma-mixture model had shown that the starting value of *b* was not critical, as long as *b* was constrained to values ≥0. In all fits, therefore, *b* was constrained to values 0 ≤ *b* ≤1, and *t*_{D}, 1/λ_{R}, or *t*_{RP} to values ≥0. In renewed fits, other starting values, particularly for the refractory parameters, were selected based on knowledge from the initial fits of the full and reduced versions of the models. The procedure was repeated iteratively and the best solution was taken.

Evaluations of the goodness-of-fit and of the randomness-of-fit of the two models were based on the vertical and horizontal differences between the sample CDF and the model CDF (the presumed population CDF). For a first analysis, we derived for each subsample a deviation measure *D* between data and fit. This measure is defined as the sum of the squared deviations, Δ^{2}, between data and fit, multiplied by the number of ISIs in the sample, *N*_{Total}, and divided by the difference between the number of ISIs in the fit, *N*_{Fit}, and the number of free parameters, *df*, in the models. Hence, *D* = Σ (Δ^{2}) × *N*_{Total}/(*N*_{Fit} − *df*). The lower this deviation measure, the better is the fit. This measure was averaged across all 399 subsamples for a given model version and for a given cut-off of short ISIs (see results). Each individual subsample is subject to noise and randomness which reduces the relative contribution of any systematic differences between a given model and the data. To emphasize the systematic differences and reduce the effects of noise and randomness in individual subsamples, we averaged, in a second analysis, the vertical (or weighted horizontal) differences between sample and model CDF at corresponding probability values. These were calculated for 5,000 probability values equally spaced between 0 and 1 by linear interpolation between the observed differences as a function of the empirical cumulative probability for each subsample. In this way, the differences are obtained at equal probabilities in different subsamples of spontaneous activity. Each difference was weighted by the square root of the number of ISIs in the sample to normalize for the fact that the magnitude of the maximum difference between a sample and a model CDF is nearly proportional to the inverse of the square root of that number (Sachs 2002). These differences were averaged across subsamples at identical probability values and divided by the mean of the square roots of the number of ISIs in the subsample. These mean residuals and their associated SE were plotted as functions of the empirical cumulative probability for both models.

##### MODEL FITTING TO A MEAN CDF.

It would be desirable to combine in a meaningful way all ISIs from all 399 subsamples into a single empirical CDF to enable the demonstration of even subtle differences between data and models. Such a procedure could be equivalent to averaging and would reduce the unsystematic noise present in individual subsamples. To obtain a mean CDF, it was necessary to transform the *x*-axis of individual CDFs before combining or averaging them. We proceeded as follows. For each of the 399 subsamples, we first calculated 4,000 ISIs for cumulative probability values equally spaced between 0 and 1 by linear interpolation between, or extrapolation from, the observed cumulative probabilities and ISIs. Second, from the equations capturing the two models (*Eqs. 3* and *8*), it can be derived that, if individual subsamples of spontaneous activity differed only in the scaling parameter λ_{E}, but not in the shape parameter *b*, the individual CDFs would simply be scaled along the (ISI − *t*_{RP})—axis and the SD of the logarithms of these differences would be independent of the value of this axis. Therefore we first subtracted from each of the 4,000 ISIs an estimate of the refractory period, *t*_{RP}, taken as 95% of the shortest ISI in the sample. We calculated the means and the SD of the natural logarithms of these differences across the 399 subsamples at corresponding cumulative probabilities. Finally, the arithmetic mean of the estimates of the *t*_{RP} across the 399 subsamples was added back to the geometric means of (ISI − *t*_{RP}), and the cumulative probability was plotted against this sum. This function is the mean CDF. It is based on 1.186 million ISIs (see results) and is thus very smooth. We fitted both models to the mean CDF, as described above.

##### SIMULATIONS.

The reliability of the fitting procedure and the parameter estimates obtained were gauged using extensive simulations, as described in the results. For illustrations and comparisons with the real data, a set of 399 simulated samples was generated according to the model that provided the better fit to the data (the gamma-mixture model). The numbers of ISIs and the mean ISIs of these simulated data sets were identical to those obtained from the real data after they had passed the stationarity analysis.

## RESULTS

### Database and generation of subsamples

Spontaneous activity was recorded from 137 AN fibers from seven owls. The durations over which this activity was recorded from individual fibers ranged from 25 s to 6 min, with a mean of 2.5 min. The total sampling time was nearly 6 h. The spontaneous discharge rates (SRs), computed from the complete recording times, ranged from about 24 to 146 spikes/s, with a mean of 72 spikes/s. They yielded ISI counts that ranged from ∼1,000 to 38,000, with a mean of 10,600, and a total count of ∼1.45 million ISIs. From these data, we generated 399 subsamples from 136 fibers, with numbers of ISIs ranging from 445 to 4,000 (mean, 2,973; total of 1.186 million ISIs), which fulfilled our strict stable-rate criteria (see methods).

The characteristic frequencies (CFs) of these fibers ranged from ∼0.5 to 9 kHz. Figure 1 shows that the SRs were inversely related to CF, a relationship that could be described by the exponential function (straight line in Fig. 1): SR = 124 × e^{−0.124 × CF}, where SR is in spikes per second and CF is in kilohertz (*r* = 0.778; *P* < 0.0001). The unexplained variance seems to be largely because of a spread of SRs at any given CF. A similar spread was seen in individual owls (different symbols in Fig. 1), suggesting that the vertical scatter is not caused by pooling of data from different owls. At any given CF, the distribution of SRs was unimodal and rather narrow, with the highest SRs being only about twice as high as the lowest SRs.

### Recording artifact affects short ISIs

When all ISIs from a given sample were sorted by their magnitudes and plotted as a CDF, a conspicuous vertical, or nearly vertical, rise of the CDF at the shortest ISIs was observed in a significant fraction of data sets. Figure 2*A* shows a representative example, with the step-like rise of the CDF highlighted in the blow-up (*inset*). Such a step results in what Köppl (1997a) has termed an “early peak” of the ISI histogram, because the bin(s) comprising that step will have by far the largest number of entries. The initial vertical rise of the CDF also translates into an early peak in the hazard function (Fig. 2*B*). Several other studies have also reported the occasional occurrence of an initial sharp peak in the hazard function, or the ISI histogram, but its origin has been rather enigmatic (Gaumond et al. 1982; Gray 1967; Klinke et al. 1994; Li and Young 1993; Miller and Wang 1993; Prijs et al. 1993; Young and Barta 1986). We have recently provided strong evidence that this peak is an artifact resulting from the common technique of recording spike times as the crossing times of a fixed trigger threshold by the electrode signal (Heil et al. 2007). As originally shown by Gaumond et al. (1982), at short ISIs, the second spike of a pair of spikes is likely to have a smaller peak amplitude than its predecessor, either as a consequence of the refractory properties of the fiber or because of the superposition of the second spike with a “hyperpolarizing” afterpotential of the first spike (see Fig. 14 in Gaumond et al. 1982; see also Fig. 5 in Siegel 1992 and Fig. 9 in Shepherd et al. 2004). The spikes recorded here also showed pronounced afterpotentials that could last for a couple of milliseconds or more (Fig. 2*B*, *inset*). The second spike with the smaller peak amplitude will therefore cross a fixed trigger threshold at a later phase than the first one, and the interval between the stored times of the two consecutive spikes will be longer than their true interval, i.e., that between corresponding phases. Gaumond et al. (1982) have estimated that these effects would lead to an apparent delay of the second spike by ∼100 μs. Consequently, the observed probability of very short ISIs is less than the true probability and that of slightly longer ISIs is higher than the true probability.

Independent evidence for the presence of such a recording artifact comes from AN fibers that discharged spikes with so-called “preferred intervals” (Köppl 1997a), even in the absence of any intentional acoustic stimulation. These preferred intervals are most likely caused by spontaneous membrane-potential oscillations associated with an electrical tuning mechanism in lower-frequency hair cells (Köppl 1997a; but see Klinke et al. 1994). An example is shown in Fig. 2, *C–F*. The preferred intervals are reflected in fluctuations of the CDF (Fig. 2*C*), but even more clearly in oscillations of the hazard function (Fig. 2, *D–F*). Generally, the hazard function describes the tendency for some event to occur at time *t* (*t* > 0, given that the last event occurred at *t* = 0) and is the probability density at time *t* renormalized by the probability that the event failed to occur before *t* (Gray 1967). In cases with preferred intervals, the hazard rate oscillates. In the particular case shown in Fig. 2, *D–F*, this occurs up to ISIs of ≥30 ms. Beyond this range, oscillations were hard to discern, possibly because of the low number of such long ISIs. A least-squares fit of a simple sine wave to the logarithm of the hazard rate between 5 and 20 ms yields an oscillating frequency of 906 Hz, close to the peak frequency of 910 Hz of a fast Fourier transform of the ISI histogram and close to the CF of 863 Hz determined from the AN fiber's tuning curve to pure tones. The gray sinusoidal lines in Fig. 2, *D–F*, show this fit, as well as its backward extrapolation to ISIs <5 ms. Enlargements of the hazard function and the sine wave fit for ISI ranges from 0 to 3 ms and from 7 to 10 ms are provided in Fig. 2, *E* and *F*. Inspection of Fig. 2*E* clearly shows that the first positive peak of the HF (at ∼1.2 ms; Fig. 2*E*, dashed vertical line) and the first negative peak (at ∼1.7 ms) are delayed, by some 250 μs, relative to the times of the positive (Fig. 2*E*, continuous vertical line) and the negative peaks expected from the backward extrapolation of the excellent fit to the HF between 5 and 20 ms (cf. Fig. 2*F*; *r* = 0.867; *n* = 7,108).

Such observations suggest that the recorded short ISIs, ≤2 ms or even more, may be subject to the recording artifact described above, even in cases where the CDF rises rather smoothly. In a related analysis, Gaumond et al. (1982) examined the responses of cat AN fibers to low-frequency tones and observed apparent phase shifts of the hazard function for ISIs even ≤3 or 4 ms. Thus to avoid bias on the model fits and their parameter estimates, we excluded, for most of the analyses described below and unless otherwise stated, all ISIs shorter than 4 ms from the fits when probing the models described in the next section. However, we also provide some results obtained with other more or less conservative cut-offs. As expected, the fraction of ISIs <4 ms increased monotonically with decreasing mean ISI (i.e., with increasing spontaneous rate), from a minimum of 0.082 for a mean ISI of 39.3 ms (25.5 spikes/s) to a maximum of 0.443 for a mean ISI of 6.59 ms (152 spikes/s). The mean fraction of ISIs <4 ms across the 399 subsamples was 0.24.

### Common framework of the models

The complete versions of both models are based on the following assumptions, in line with nearly all models of AN fiber activity. *1*) Excitation to an AN fiber is provided by a stochastic point process. Hence, excitatory events occur at random times. *2*) Each excitatory event results in a spike unless the fiber is in a refractory state. *3*) The refractory state is composed of an absolute refractory period of duration *t*_{D}, followed by a relative refractory period during which the probability *r*(*t*) of a spike occurring given an excitatory event increases exponentially from 0 to 1 with time since the previous spike, with a time-constant 1/λ_{R} (1) where *t* is the time since the previous spike. This refractory function has been implemented in many previous modeling studies (Heil et al. 2007; Li and Young 1993; Lütkenhöner et al. 1980; Meddis 2006; Meddis and O'Mard 2005; Prijs et al. 1993; Schoonhoven et al. 1997; Schroeder and Hall 1974; Young and Barta 1986), and its shape is in excellent agreement with physiological data (Brown 1994; Miller et al. 2001; Morsnowski et al. 2006). The fiber enters a refractory state after each spike and the refractory memory extends back in time only to the last spike. *4*) Excitatory events and refractoriness are independent, i.e., the fiber's refractory state has no influence on the generation of excitatory events.

The most common model of this type used to account for ISI distributions during spontaneous activity contains a homogeneous Poisson process for excitation (Bi 1989; Carney 1993; Geisler 1981, 1998; Geisler et al. 1985; Javel 1986; Kiang et al. 1965; Krishna 2002; Kuhlmann et al. 2002; Li and Young 1993; Lütkenhöner et al. 1980; Manley and Robertson 1976; Miller and Wang 1993; Molnar and Pfeiffer 1968; Schmich and Miller 1997; Schroeder and Hall 1974; Young and Barta 1986; Zhang et al. 2001; for reviews, see Delgutte 1996; Mountain and Hubbard 1996). We have shown, however, that in the cat, this cannot be the case, particularly in combination with physiologically plausible refractoriness (Heil et al. 2007). Instead, the ISI distributions were both better and more accurately described by nearly 2 orders of magnitude with a homogeneous stochastic process of excitation in which the distribution of intervals between excitatory events is a mixture of an exponential and a gamma distribution with shape factor 2, both with the same scaling parameter. Notably, this solution emerged from an attempt to rescue the Poisson model by modifying the refractory function (Heil et al. 2007).

In the following, we therefore examine the standard (Poisson) model and the gamma-mixture model that was adequate for the cat. To keep both models simple, the oscillations of the hazard rate around a mean rate (as in Fig. 2, *D–F*) are ignored. Modeling their complex properties is beyond the purpose of this study and did not affect the conclusions of this paper.

### Poisson model: excitatory events are exponentially distributed

According to this Poisson model, the CDF of the intervals between two consecutive excitatory events [the interevent interval (IEI)] is given by the exponential distribution (2) where *t* (in s) denotes the duration of the IEIs and λ_{E} (in s^{−1}) is the scaling factor or the rate of the homogeneous Poisson process. Figure 3 shows the PDF and the HF of the IEIs of this process. The PDF is a straight line of negative slope on a logarithmic probability density axis, and the hazard rate is constant and equal to λ_{E}.

In this scenario, the ISIs should be distributed according to a general-gamma distribution (Heil et al. 2007; Lehmann 2002; Li and Young 1993; Young and Barta 1986), whose CDF is given by (3)

For ISIs much longer than the sum of the absolute refractory period, *t*_{D}, and of the mean duration of the relative refractory period, 1/λ_{R}, i.e., (*t* − *t*_{D}) ≫ 1/λ_{R}, and provided the mean duration of relative refractory period is short compared with the mean interval between excitatory events, i.e., 1/λ_{R} ≪ 1/λ_{E}, *Eq. 3* is approximated by the exponential distribution (4) where *t*_{RP} = *t*_{D} +1/λ_{E} × ln[λ_{R}/(λ_{R} − λ_{E})] ≅ *t*_{D} +1/λ_{R} (see appendix A for a formal derivation of this approximation). Hence, parameter *t*_{RP} is a close approximation of the sum of the absolute refractory period, *t*_{D}, and of the mean duration of the relative refractory period, 1/λ_{R}. This sum will henceforth be referred to as the mean refractory period. *Eq. 4*, which has only two free parameters, can also be conceived of as a Poisson model with only a (prolonged) absolute refractory period. It will be referred to as the reduced version of the Poisson model. The mean μ_{ISI} of an ISI distribution of the Poisson model is given by (5)

### Gamma-mixture model: the distribution of IEIs is a mixture of an exponential and a gamma distribution with shape factor 2

This model assumes that the distribution of IEIs is a mixture of an exponential distribution (which is a gamma distribution with shape factor β = 1) and a gamma distribution with shape factor β = 2, where both components have the same scaling factor, λ_{E}. Their relative proportions in the mixture will be denoted as (1 − *b*) and *b*, respectively. The CDF of the IEIs is thus given by (6)

This mixture can emerge from a primary homogeneous Poisson point process with rate λ_{E}, if some of these primary events are removed by a subsequent process or fail to trigger secondary events further down a cascade of precursor steps that might be required to produce the ultimate excitatory events, i.e., those events which trigger spikes unless they fall into the refractory period. If every other event generated by the primary Poisson process is removed, the intervals between the remaining events have a gamma distribution with shape factor β = 2 and scaling factor λ_{E} (Cox 1962; Koch 1999; Stein 1965). If this happens only intermittently, the distribution of the intervals between the remaining events will be a mixture of a gamma distribution with β = 2 and an exponential distribution of common scaling factor. In other words, primary events can be removed from the Poisson process, if the following simple rule is obeyed: never remove two or more consecutive events.

Figure 3 (gamma-mixture model) shows the PDFs and the HFs of a gamma distribution with shape factor β = 2 and scaling factor λ_{E}, of an exponential distribution (a gamma distribution with shape factor β = 1) with the same scaling factor, and of a mixture of the two with fractions of *b* and (1 − *b*), respectively. The PDF of the gamma distribution with β = 2 is a nonmonotonic function that reaches its peak at 1/λ_{E}. If the probability density is plotted on a logarithmic axis, the steepness of the descending limb of the PDF increases as *t* increases. The PDF for the mixture of distributions lies somewhere in between those for the gamma and the exponential distribution, with the relative distances to the latter two depending on *b*. This also holds for the HF. The HFs rise monotonically from 0 for the gamma distribution and from (1 − *b*) λ_{E} for the mixture to an asymptotic limit of λ_{E}.

Assuming the same refractory function (*Eq. 1*) as in the Poisson model, the CDF of the gamma-mixture model is given by Heil et al. (2007) (7)

The first summand, with the weighting factor of (1 − *b*), describes the distribution of that fraction of ISIs that results from the exponentially distributed intervals between excitatory events (cf. *Eq. 3*), and the second summand, with the weighting factor of *b*, the distribution of the other fraction of ISIs that results from the gamma distributed intervals between excitatory events.

As with the Poisson model, the complexity of the gamma-mixture model can be reduced, particularly when short ISIs are excluded from the fits, a procedure necessitated by the recording artifact. It can be shown (see appendix B) that for ISIs much longer than the mean refractory period, i.e., (*t* − *t*_{D}) ≫ 1/λ_{R}, and provided the mean duration of the relative refractory period is short compared with the mean interval between excitatory events, i.e., 1/λ_{R} ≪ 1/λ_{E}, *Eq. 7* is approximated by the following simpler equation (8) where *t*_{RP} ≅ *t*_{D} +1/λ_{R} as before. *Equation 8* can also be conceived of as an independent model, viz., a gamma-mixture model with only a (prolonged) absolute refractory period. This reduced version of the gamma-mixture model has only three free parameters.

The mean μ_{ISI} of an ISI distribution of the gamma-mixture model is given by (9) For *b* = 0, the gamma-mixture model is identical with the Poisson model.

According to the gamma-mixture model, the modest slow monotonic rise of the hazard function over at least several tens of milliseconds (Fig. 2*B*), a phenomenon even more clearly seen in mammals (Gaumond et al. 1982, 1983; Gray 1967; Heil et al. 2007; Li and Young 1993; Prijs et al. 1993), originates from the fraction of gamma-distributed IEIs of the stochastic process providing the excitation (Fig. 3). The refractory properties of the AN fiber, if brief, only contribute to shaping the initial rapid rise of the HF. With the assumption of a Poisson process providing the excitation, the slow rise of the HF would have to be attributed to the refractory properties of AN fibers.

### Reduced gamma-mixture model is superior to the full Poisson model

The reduced version of the gamma-mixture model (*Eq. 8*) can be viewed as an independent model, viz., a gamma-mixture model with only a (prolonged) absolute refractory period. It has the same number of free parameters (3) as the full version of the Poisson model (*Eq. 3*). Therefore these two versions of the two different models can be compared directly when examining their success of capturing the empirical ISI distributions. We fitted these model versions to all 399 subsamples, either including all ISIs in the fit or excluding ISIs shorter than some specified value. Excluding short ISIs diminished the influence of the recording artifact on the fits. For each subsample, we first derived a deviation measure *D* (see methods). The lower this deviation measure, the better is the fit. Figure 4 *A* plots the average of these measures across the 399 subsamples as a function of the lower bound of ISIs included in the fit. The figure clearly shows that the reduced version of the gamma-mixture model (•) is considerably better than the full version of the Poisson model (○), irrespective of whether the entire sample CDFs were fitted or whether short ISIs were excluded. The difference between the means of the deviation measures for the two models increases when the short ISIs are excluded, ≤3 or 4 ms, presumably because the influence of the recording artifact decreases. For all lower bounds, the gamma-mixture model is significantly better than the Poisson model (Wilcoxon tests for matched pairs yielded *P* < 10^{−12} for all lower bounds examined). The distribution of the 399 ratios of the deviation measures obtained with the Poisson and the gamma-mixture models for the lower bound of 4 ms is shown in Fig. 4*B*. Although a number of ratios are very close to 1, the distribution is smooth and provides no grounds for the assumption of separate fiber populations, with some fraction being better described by the gamma-mixture model and the other fraction being better described by the Poisson model. The simpler assumption that the fibers rather form a homogeneous population is further supported by simulations. The gray line shows the corresponding distribution of 399 ratios obtained from a set of simulated data, generated according to the reduced version of the gamma-mixture model with a fixed refractory period *t*_{RP} of 0.403 ms and a fixed *b* of 0.225. The two ratio distributions are very similar. The maximum vertical differences between the distributions of values from real and simulated data (0.0872) is smaller than the critical difference *D*_{α} of 0.0963 for α = 0.05. The critical difference is given by *D*_{α} = *K*_{(}_{α}_{)} × [(*n*_{1} + *n*_{2})/(*n*_{1} × *n*_{2})]^{0.5}, where *K*_{(}_{α}_{)} is a constant that only depends on the significance level. For α = 0.05, *K*_{(}_{α}_{)} = 1.36 (Kolmogorov-Smirnov tests; Sachs 2002). Here, *n*_{1} = *n*_{2} = 399. This means that the null hypothesis of equality of the two distributions (from real and from simulated data) of ratios of deviation measures should not be rejected. In summary, these data show that the gamma-mixture model is superior to the Poisson model in its capability to describe the data. It seems that the differences in model structure, and not model complexity, account for these different capabilities. However, our finding that the reduced version of the gamma-mixture model provides a better fit than the full version of the Poisson model can also be viewed in the following way. It could be interpreted as showing that a “remove events, but never two in a row, from a Poisson process” procedure is better in explaining the ISI distributions than a “remove events from a Poisson process via a relative refractory period” procedure. Finally, the ISI distributions of all AN fibers, and not just some fraction, are better described by the gamma-mixture model.

### Reduced version of the gamma-mixture model can capture the empirical coefficients of variation and the skew, whereas the full version of the Poisson model cannot

For each subsample and with all ISIs included, we computed the CV (defined as the ratio of the SD to the mean ISI) and the skew. Figure 5 plots the mean (intersection of thick bars, located at a CV of 0.919 and a skew of 1.775) and 4 SD of the CV and of the skew across the 399 subsamples (thick bars). To examine whether these empirical measures can be reproduced by the models, we generated numerous sets of 399 finite samples from theoretical CDFs calculated according to the Poisson model and the gamma-mixture model, using systematic variations of the model parameters. The numbers of ISIs for each of the 399 such samples were chosen such that they were identical to their 399 counterparts of real subsamples. Furthermore, the values of λ_{E} were chosen, according to *Eqs. 5* and *9*, such that the resulting mean ISIs were identical to those of their 399 real counterparts. For a given set of 399 generated samples, the values of *t*_{D} (or *t*_{RP}), 1/λ_{R}, and *b* were held constant, but were varied systematically between sets. For each set, the resulting CV and skew values were averaged across the 399 samples. Informative data are shown in Fig. 5. The gray solid line with open circles labels the trajectory in the CV-skew plane that is obtained with the Poisson model (*b* = 0), with *t*_{D} = 0 ms and with 1/λ_{R} increasing from 0 (top right open circle) to 3 ms (bottom left open circle; the other circles are placed at intervals of 1/λ_{R} of 0.5 ms). This line marks the upper boundary of CV-skew combinations that can be achieved with physiologically possible realizations of this model, namely when the absolute refractory period is equal to or greater than zero. Clearly, the empirical combination of CV and skew falls above this upper limit. The empirical combination can be achieved with the Poisson model only when the absolute refractory period parameter of the interval distribution is set to the negative value of −1.394 ms (gray dashed line). Hence, the empirical combination of CV and skew, even when considering the confidence intervals (CIs) (see thick bars in Fig. 4), cannot be reproduced by the Poisson model with physiologically possible parameter values. The null hypothesis that the 399 data samples are drawn from a population generated according to the Poisson model has to be rejected. The gamma-mixture model, on the other hand, can easily reproduce the empirical combination with physiologically possible parameter values. The black solid line with open triangles marks the trajectory of the gamma-mixture model with *t*_{D} = 0 ms and also with 1/λ_{R} = 0, and with *b* increasing from 0 (top right triangle, identical with gray open circle) to 0.30 (bottom left triangle; the other triangles mark intervals of *b* of 0.05). This line marks the upper boundary of CV-skew combinations that can be achieved with physiological realizations of the gamma-mixture model, namely when the refractory period is equal to or greater than zero. The empirical CV-skew combination falls below this upper limit. This combination can be reproduced when a physiologically plausible refractory period is included in the model. The black line with filled triangles (marking intervals of *b* of 0.05), which is parallel to that with open triangles, represents the trajectory for an absolute refractory period *t*_{D} of 0.543 ms (note that here *t*_{D} = *t*_{RP} because a relative refractory period is not included in this model; cf. *Eq. 8*). It intersects the empirical CV-skew combination at a value of *b* = 0.237. Finally, the vertical black line with oblique crosses shows the trajectory for this model with *b* = 0.237 and with *t*_{D} increasing from 0 (top cross) to 2 ms (bottom cross; other crosses mark intervals of *t*_{D} of 0.5 ms). It intersects the empirical CV-skew combination at *t*_{D} = 0.543 ms.

In summary, these data underscore the superiority of the gamma-mixture model over the Poisson model. The former easily accounts for the combination of CV and skew seen in the data, but the Poisson model fails to do so with physiologically possible parameter values. Furthermore, the analysis of the average CV and skew of the real and simulated ISI distributions performed here constitutes one method by which the values of the model parameters *b* and of *t*_{RP} can be estimated. The average empirical CV-skew combination is reproduced with *b* = 0.237 and *t*_{RP} = 0.543 ms.

### Validity of the approximations

For a more detailed analysis, we excluded ISIs <4 ms from the fits for the reasons given above, and in the following compare the results obtained with the full version of the Poisson model (*Eq. 3*) to those obtained with the reduced version of the gamma-mixture model (*Eq. 8*). To show the validity of the approximations of the full by the reduced models under these conditions, we fitted each of the 399 subsamples of spontaneous activity with the full and with the reduced versions of both models for ISIs >4 ms. Figure 6 provides comparisons of the mean and the associated SD of the deviation measures (Fig. 6*A*) and of the parameter estimates (Fig. 6, *B–D*) obtained from the fits of both versions of both models to each of the 399 subsamples. The figure allows appreciation of the fact that deviation measures and parameter estimates were very similar for the reduced and full versions of each model, respectively. For the reduced and full versions of the Poisson model (Fig. 6, 2 left bars), the average deviation measures and associated SD were 0.0798 ± 0.0039 and 0.0790 ± 0.0038, respectively. The average and SD of the mean interval between excitatory events, 1/λ_{E}, were 13.791 ± 0.246 and 13.784 ± 0.246 ms for the reduced and full versions, respectively, a mean difference of only 7 μs. The average and SD of the mean refractory periods, i.e., *t*_{RP} and (*t*_{D} + 1/λ_{R}), were 0.646 ± 0.016 and 0.637 ± 0.016 ms for the reduced and full versions, respectively, a mean difference of 9 μs.

Similar results were obtained when comparing the reduced and the full versions of the gamma-mixture model (Fig. 6, 2 right bars in each panel). The means and SD of the deviation measures were 0.0468 ± 0.0023 for both the reduced and the full versions. The average and SD of 1/λ_{E} were 11.556 ± 0.193 and 11.560 ± 0.193 ms, a mean difference of 4 μs. The average and SD of the mean refractory periods were 0.363 ± 0.016 and 0.361 ± 0.016 ms, a mean difference of 2 μs. The average and SD of *b* were 0.202 ± 0.007 and 0.201 ± 0.007 for the reduced and the full versions, respectively.

In summary, these analyses clearly show the validity of the approximations of the full models by the reduced models. The analyses therefore also show that a relative refractory period of barn owl AN fibers must be rather short, because it has no appreciable influence on the distribution of ISIs >4 ms. Furthermore, they show that the lack of parameter *b* in the full version of the Poisson model cannot be compensated for by the relative refractory period: the reduced version of the gamma-mixture model is superior (compare the heights of the 2 center bars in Fig. 6*A*; see also Fig. 4). The nature of the differences in the fits of these two models is explored in more detail in the next sections.

### Comparing the Poisson and the gamma-mixture models

##### A REPRESENTATIVE EXAMPLE.

Results obtained with fits of the full version of the Poisson model (3 free parameters) and of the reduced version of the gamma-mixture model (also 3 free parameters) to the ISI distribution, excluding ISIs <4 ms, of one representative example are shown in Fig. 7. Figure 7*A* shows the empirical CDF along with the fit of the full version of the Poisson model (*Eq. 3*; gray line). It seems to provide a good fit to the data, except at ISIs shorter than ∼3 ms. Here, the Poisson model underestimates the observed probabilities. If the Poisson model were correct, but short ISIs were prolonged because of the recording artifact explained above, the model fit should have overestimated, and not underestimated, the observed probabilities of the short ISIs. In this example, 1/λ_{R} is close to zero, whereas *t*_{D} is long compared with 1/λ_{R}, but this was not true in general. With ISIs <4 ms excluded, estimates of 1/λ_{R} ranged from 0 to ≤1.53 ms (mean, 0.49 ms), and those of *t*_{D} from 0 to 1.37 ms (mean, 0.15 ms). Figure 7*B* shows the corresponding residual probabilities (data minus fit) plotted as a function of ISI. The residuals show that the Poisson model, despite its seemingly good fit for ISIs >4 ms (Fig. 7*A*), overestimates the probabilities of medium ISIs (between ∼6 and 12 ms) and underestimates those of long ISIs (>12 ms). Across all ISIs, the residuals form an M-shaped pattern. Figure 7*C* shows the residuals in the time domain after weighting with the survival function. The square of this measure was one of the factors in the cost function minimized during initial steps of the fitting routine (see methods). As expected, these time residuals show the inverse pattern of the probability residuals (cf. Fig. 7*B*): the short ISIs are shorter than predicted by the model, intermediate ISIs longer, and long ISIs again shorter than predicted.

Figure 7, *D–F*, shows the analogous results obtained with the fit of the reduced version of the gamma-mixture model (*Eq. 8*), which for this example yielded an estimate of the presumed fraction of gamma-distributed IEIs of *b* = 0.346. The fit at ISIs >4 ms seems very good (Fig. 7*D*). The model overestimates the probabilities of very short ISIs, a discrepancy expected if the model were correct and short ISIs were prolonged because of the recording artifact. Inspection of the residuals in the vertical (Fig. 7*E*) and horizontal domain (Fig. 7*F*) shows that, over the range of ISIs fitted, the residuals are considerably smaller than with the Poisson model (cf. Fig. 7*E* with *B* and *F* with *C*).

##### GROUP DATA, GOODNESS-OF-FIT, AND RANDOMNESS OF FIT.

The general validity of the results just described for the representative example is emphasized by the mean vertical and mean horizontal differences across all 399 subsamples (see methods). This averaging reduces the noise and randomness in the individual subsamples and thus emphasizes systematic differences between model and sample CDFs. These are more readily suited for a comparison of the two models and for an evaluation of their randomness-of-fit than individual data sets. For each model, the individual parameter estimates that yielded the best fits to a given subsample were used to generate the mean differences. Figure 8 (*top 4 panels*) plots the mean vertical and the mean weighted horizontal differences between sample CDFs and model CDFs as functions of the empirical cumulative probability, for the Poisson model (Fig. 8, *A* and *B*) and the gamma-mixture model (Fig. 8, *E* and *F*).

For the Poisson model, the mean vertical and horizontal differences show a pronounced M-shaped or W-shaped pattern, respectively, and the mean differences substantially exceed the confidence intervals (Fig. 8, gray lines). The CIs start to widen with decreasing probability because ISIs <4 ms were excluded from the fits. For the gamma-mixture model, the differences are considerably smaller. Only at low probabilities (<0.3; corresponding to short ISIs) are the observed mean probabilities drastically lower, or, equivalently, the observed ISIs longer, than predicted by the model. As mentioned above, this mismatch is expected if ISIs were prolonged by the recording artifact. Nevertheless, even at higher probabilities, the mean vertical and horizontal differences for the gamma-mixture model do exceed the confidence limits. We examined the possibility that these residuals that remain unexplained by the gamma-mixture model might be caused by AN fibers showing preferred intervals. However, this was not the case. The patterns of residuals were very similar for the groups of subsamples from fibers with (*n* = 91) and without preferred intervals (*n* = 308) and were both similar to those shown in Fig. 8.

##### QUANTITATIVE COMPARISONS OF THE GOODNESS-OF-FIT OF THE TWO MODELS.

For quantitative comparisons of the goodness-of-fit of the two models, we restricted ourselves to measures taken from those curves for probabilities >0.443. This is done because the cumulative probability reached at a given ISI varies with the mean spontaneous rate; the higher the rate, the higher the probability. At an ISI of 4 ms, corresponding to the lower limit used for inclusion of data points in each fit, the highest cumulative probability reached by any subsample was 0.443. For probabilities >0.443, the maximum absolute vertical difference was 2.9 times larger, the maximum horizontal weighted difference 2.2 times larger, the average of the squared mean vertical differences 10.0 times larger, and the average of the squared mean horizontal differences 11.8 times larger for the Poisson model than for the gamma-mixture model. These comparisons show, even more clearly than Figs. 4 and 6, that the reduced version of the gamma-mixture model provides a much better fit to the data than the full version of the Poisson model.

### Estimates of parameters returned by the gamma-mixture model

Because the gamma-mixture model provides the better fit to the empirical CDFs, it is interesting to examine the estimates of the three parameters returned by the fits. We examine those derived from fits to ISIs >4 ms. As explained above, parameter *t*_{RP} is a close approximation of the mean refractory period, i.e., of the sum of the absolute and of the mean duration of the relative refractory periods. Estimates of *t*_{RP}, which were constrained to values ≥0, varied between 0 and 1.55 ms, with a mean and SD of 0.363 ± 0.317 ms across the 399 subsamples. Values of *t*_{RP} very close to 0 were not uncommon and were encountered across a wide range of mean ISIs. These physiologically implausible values are not associated with a special subgroup of AN fibers but are likely the result of the chance distribution of ISIs in a given subsample. This conclusion is supported by further analyses that showed that estimates of *t*_{RP} taken from different subsamples of spontaneous activity of a given individual AN fiber could cover a wide range including values close to zero. It is also supported by the results obtained from simulations. The estimates of *t*_{RP} were independent of the mean ISI (*r* = −0.0542; *P* = 0.28053; *n* = 399) and independent of CF (*r* = −0.0487; *P* = 0.33829). When the estimates of *t*_{RP} from different subsamples of the same fiber, which are not independent, were first averaged before calculating the mean and SD across the population of AN fibers, rather than across that of the subsamples, a similar mean and SD of *t*_{RP}, viz. 0.365 ± 0.266 ms, were obtained. These estimates, which are plotted against the corresponding mean ISI in Fig. 9*A*, were also independent of the mean ISI (*r* = −0.0979; *P* = 0.25693; *n* = 136) and independent of CF (*r* = −0.0444; *P* = 0.61725). The estimates of *t*_{RP} are subject to a modest bias (see *Simulations* below), so their mean of 0.365 ms should not be taken as the best estimate of the mean refractory period.

Estimates of parameter *b*, which were also constrained to values ≥0, varied between 0 and ∼0.52, with a mean and SD of 0.202 ± 0.140. Values of *b* very close to 0 were also encountered across a wide range of mean ISIs. Again, these values are unlikely associated with a special subgroup of AN fibers but are more likely the result of the chance distribution of ISIs in a given subsample. This conclusion is supported by the finding that estimates of *b* taken from different subsamples of spontaneous activity of a given individual AN fiber could cover a wide range including values close to zero and by the results obtained from simulations. Nevertheless, *b* and mean ISI were positively correlated (*r* = 0.2504; *P* = 4.06 × 10^{−7}; *n* = 399), indicating that *b* tends to decrease as spontaneous rate increases. Because spontaneous rate increases with decreasing CF (Fig. 1), *b* and CF were also positively correlated (*r* = 0.1607; *P* = 0.00149; *n* = 388). When the dependent estimates of *b* from different subsamples of the same fiber were first averaged, a mean and SD of *b* of 0.211 ± 0.116 across the population of AN fibers was obtained. These estimates, which are plotted against the corresponding mean ISI in Fig. 9*B*, were also positively correlated with mean ISI (*r* = 0.3339; *P* = 7.11 × 10^{−5}; *n* = 136; Fig. 9*B*) and with CF (*r* = 0.2555; *P* = 0.00347).

Figure 9*C* shows a plot of the time constant of the excitatory process, 1/λ_{E}, against the mean ISI. Rearranging *Eq. 9* shows that 1/λ_{E} should be related to the mean ISI by (10) If *b* and *t*_{RP} were identical across samples, the data points ought to fall on a straight line with a slope of 1/(1 + *b*) and an intercept of −*t*_{RP}/(1 + *b*). As expected from the scatter of the estimates of *t*_{RP} and *b* and the correlation of the latter with mean ISI, the data points scatter around, and fan out from, the line computed with the means of *b* and *t*_{RP} (Fig. 9*C*, dashed line).

### Simulations

To further explore the reasons behind the observation that the mean residuals of the gamma-mixture model exceed the confidence limits (Fig. 8, *E* and *F*) and to gauge the reliability of the parameter estimates, we performed extensive simulations of ISI distributions generated according to the simplified version of that model (*Eq. 8*). Each ISI was generated as the sum of a refractory period, *t*_{RP}, and one interval (for the fraction of exponentially distributed ISIs) or two intervals (for the fraction of gamma distributed ISIs) drawn from an exponential distribution with mean 1/λ_{E} (cf. *Eq. 9*). The exponentially distributed intervals were computed from the uniformly distributed random number *z*, where 0≤ *z* ≤1, as −1/λ_{E} ln(1 − *z*). In a first round of simulations, we generated ISI distributions, keeping *t*_{RP} and *b* fixed at the means derived from the fits of the gamma-mixture model to the real data, viz. *t*_{RP} at 0.363 ms and *b* at 0.202 for all simulations. We simulated each of the 399 real subsamples of spontaneous activity 15 times using that value for λ_{E} (calculated by the method-of-moments from *Eq. 9*) that would be expected to yield a mean ISI identical to that of the real subsample to be simulated. Furthermore, each simulated sample contained the same number of ISIs as the real subsample to be simulated. Each of the 15 × 399 simulated samples was fitted with the Poisson model and with the gamma-mixture model, excluding all ISIs <4 ms just as for the real data. The means of the estimated parameters *t*_{RP} and *b* derived from these fits were ∼10% smaller than the values used to generate the data, reflecting a slight bias whose possible origin is pursued in the next section. In a second round of simulations, we compensated for this bias by multiplying the values with the factors necessary to obtain mean estimates of these parameters very close to those obtained from the real data. The new values used were *t*_{RP} = 0.403 ms and *b* = 0.225.

Figure 8 (*bottom 4 panels*) plots the mean vertical and the mean weighted horizontal differences between the first set of 399 simulated sample CDFs out of 15 sets from this second round of simulations and the model CDFs as functions of the empirical cumulative probability, for the full version of the Poisson model (Fig. 8, *C* and *D*) and the reduced version of the gamma-mixture model (Fig. 8, *G* and *H*). For the Poisson model, the mean vertical differences and mean horizontal differences show similar M- and W- shaped patterns, respectively, as were observed for the real data (cf. Fig. 8, *A* and *B*). Only the magnitudes of the differences at low probabilities are considerably smaller for the real data than for the simulated data. This difference is readily explained by the presence of the recording artifact in the real data and its absence in the simulated data. The artifact reduces the mismatch of the real data with the Poisson model at low probabilities (or, equivalently, at short ISIs). For the gamma-mixture model, the mean vertical and the mean horizontal differences between the simulated sample CDFs and the model CDFs fluctuate unsystematically around zero and fall nearly entirely within the confidence limits (Fig. 8, *G* and *H*). Only at very low probabilities (<0.015) are the mean vertical differences considerably more negative than the lower confidence limit (Fig. 8*G*). This is a boundary effect resulting from the restriction that probabilities cannot be smaller than zero. Hence, at these very low probabilities, there is pronounced asymmetry of possible vertical differences. There is no such asymmetry of possible horizontal differences because the refractory period shifts the ISIs far enough away from the boundary of zero, and consequently, the mean horizontal differences stay within the confidence limits even at very low probabilities (Fig. 8*H*). When the mean vertical and horizontal differences of the simulated data to the gamma-mixture model (Fig. 8, *G* and *H*) are compared with the corresponding differences of the real data (Fig. 8, *E* and *F*), it is obvious that the latter are much larger. Thus the latter cannot be attributed to the fitting routine. The large differences at low probabilities in the real data, and their absence in the simulated data (with the boundary effect on the vertical differences set aside), are consistent with the effects expected from the recording artifact.

As expected, the estimates of *t*_{RP} and of *b* obtained from the fits of the gamma-mixture model to the sets of 399 simulated sample CDFs were uncorrelated with the mean ISI. For the first simulated set, *r* = −0.039; *P* = 0.44193 for *t*_{RP} and *r* = −0.017; *P* = 0.73085 for *b*. The means and SD of these estimates were *t*_{RP} = 0.381 ± 0.262 ms and *b* = 0.201 ± 0.124. The means are thus close to the values obtained from the 399 subsamples of real data (which were 0.363 ms for *t*_{RP} and 0.202 for *b*), but the SD is smaller (which were 0.317 ms for *t*_{RP} and 0.140 for *b*). The latter result is consistent with the idea that *t*_{RP} and *b* may vary between AN fibers in the owl.

Analyses and comparisons of the distributions of the estimates from 15 sets of 399 such simulated data and the real data (Fig. 10) support this idea. The corresponding distributions are significantly different, although their total widths are roughly similar, and their means very similar (being 0.370 vs. 0.363 ms and 0.209 vs. 0.202 for *t*_{RP} and *b* from simulated and real data, respectively). For both parameters, the maximum vertical differences between the distributions of values from real and simulated data (0.0961 for *t*_{RP} and 0.0886 for *b*) were larger than the critical difference *D*_{α} of 0.0843 for α = 0.01. It is computed as *D*_{α} = *K*_{(}_{α}_{)} × [(*n*_{1} + *n*_{2})/(*n*_{1} × *n*_{2})]^{0.5}, where, for α = 0.01, *K*_{(}_{α}_{)} = 1.63 (Kolmogorov-Smirnov tests; Sachs 2002), and in our case, *n*_{1} = 399 and *n*_{2} = 15 × 399 = 5,985. This means that the null hypothesis of equality of the two distributions (of estimates from real and from simulated data) has a likelihood of being correct of *P* = 0.0019 for *t*_{RP} and *P* = 0.0055 for *b* and should therefore be rejected. As expected, estimates derived from simulated data generated with variable *t*_{RP} and variable *b* were more broadly distributed than those from data generated with fixed values of these parameters (data not shown).

Of course, we cannot exclude the possibility that the differences in the distributions of parameter estimates derived from the real and the simulated data (Fig. 10) arise because the gamma-mixture model is not a perfect description of the ISI distributions in barn owl primary afferents, not even of their ISI distributions. Nevertheless, the differences in the distributions are in line with the idea, derived from the correlation of *b* estimates with SR and with CF, that parameter *b*, but possibly also *t*_{RP}, varies between primary afferents in the barn owl.

### Possible origin of the bias

We performed further analyses on simulated data to address the possible origin of the bias. First, the fitting procedure did not always underestimate *t*_{RP} and *b*. Across the fits performed on the 5,985 simulated data sets generated with *t*_{RP} = 0.403 ms and *b* = 0.225, *t*_{RP} was underestimated in 53.72% and overestimated in 46.28% of cases, whereas *b* was underestimated in 49.96% and overestimated in 50.04% of cases (see also Fig. 10). Second, we sorted the 5,985 estimates of each parameter according to the number of ISIs in the fit and computed the mean and the SD of both parameter estimates in six nonoverlapping bins of numbers of ISIs. Each bin contained ∼1,000 values. We plotted the means and the associated SD against the mean number of ISIs in the samples in each of those six bins. For both parameters, SD decreased systematically with increasing number of ISIs, showing that the distributions of parameter estimates are narrower the larger the number of ISIs. Conversely, the corresponding means increased. They approached the parameter values used to generate the data. Hence, as the number of ISIs increases, the estimates become more reliable and the bias diminishes. This is further confirmed by the analysis of the mean CDF computed from these simulated data, and as described in the next section.

Because sample distributions approach the theoretical distribution as the number of observations in the samples increases (the Glivenko-Cantelli theorem; Sachs 2002), our fitting procedure should yield bias-free parameter estimates when the ISIs are drawn from the theoretical distribution. This was indeed the case. Fits of such distributions returned the same parameter values as those used to generate the data (within the accuracy limits of the computer), a control also performed in our previous study of cat AN data (Heil et al. 2007). Hence, our procedure can retrieve the model parameters accurately. This suggests that the small bias ultimately results from the stochasticity of the ISIs.

### Analysis of mean CDFs

As a complementary method for deriving the mean residuals (Fig. 8) and parameter estimates and for checking the conclusion that the observed differences between the residuals in the real and the simulated data with the gamma-mixture model may be caused by the recording artifact, we calculated mean CDFs from both real and simulated data. Figure 11 *A* (black line) shows the mean CDF, calculated as explained in methods, for the real data. The gray line shows the best fit of the gamma-mixture model (to ISIs >4 ms). The vertical and weighted horizontal differences of this fit were almost identical to the functions shown in Fig. 8, *E* and *F*, respectively, and are therefore not shown. This shows that both forms of averaging lead to a valid measure of the systematic differences between data and model. This conclusion is further supported by a comparison of the parameter estimates. The fit of the reduced version of the gamma-mixture model to the mean CDF (to ISIs >4 ms) yielded estimates of *t*_{RP} = 0.535 ms and of *b* = 0.224. The means of the estimates derived from the fits to individual empirical CDFs (to ISIs >4 ms) were slightly lower, viz. *t*_{RP} = 0.363 ms and *b* = 0.202, but these values correspond to values of about *t*_{RP} = 0.403 ms and *b* = 0.225, when corrected for the bias of the fitting procedure. These numbers are also similar to the estimates derived from the CV-skew analysis, which yielded *t*_{RP} = 0.543 ms and of *b* = 0.237.

The *inset* of Fig. 11*A* shows a blow-up of the initial portions of the empirical mean CDF and the fit. Because of the recording artifact, the short mean empirical ISIs are longer than predicted by extrapolation of the gamma-mixture model. The maximum difference is ∼300 μs, three times larger than the 100 μs estimated by Gaumond et al. (1982).

Figure 11*C* shows the corresponding mean CDF for simulated data, generated with *t*_{RP} = 0.403 ms and with *b* = 0.225, and the nearly perfect fit of the gamma-mixture model, which yielded estimates of *t*_{RP} = 0.414 ms and of *b* = 0.226. The differences between simulated data and fit over the extrapolated range of short ISIs are very small (see *inset*).

Figure 11*B* shows, for the real data, the SD of the logarithms of the differences between the ISIs and the *t*_{RP}, plotted as a function of ISI. As explained in methods, if the CDFs from different subsamples of spontaneous activity would differ only in their scaling parameter, this SD should be constant across all ISIs > *t*_{RP} (and 0 for ISI ≤ *t*_{RP}). However, the SD obtained from the real data rises rather steeply to reach a pronounced maximum at ∼0.9 ms and declines quickly until ∼3 ms; beyond that, the decline is shallower. A minimum of ∼0.34 is reached at ∼50 ms, beyond which the SD rises again. The mean SD (for cumulative probabilities from 0.1 to 0.9) is 0.368 (horizontal line).

Figure 11*D* (black) shows the SD data from the corresponding simulations. The function rises very steeply from almost zero at 0.403 ms (the value of *t*_{RP} used to generate the data) to quickly reach an essentially constant plateau of 0.35 maintained out to very long ISIs. Only beyond 70 ms does the SD increase slightly. The mean SD is 0.350. A very similar function was obtained from simulated data generated with a variable *b* (data not shown). When compared with the SD function obtained from the real data (cf. Fig. 11*B*), the most striking difference is the presence of the pronounced peak in the SD function from the real data and its absence from the simulated data. Given the location of this peak, with its maximum at 0.9 ms, it most likely results from the recording artifact. This artifact leads to a violation of the assumption of scaled CDFs.

## DISCUSSION

### Poisson and gamma-mixture models

In this study, we analyzed aspects of the distribution of ISIs of spontaneous activity of AN fibers in the barn owl. We modeled the distributions of ISIs as resulting from the interaction of a refractory period with a sequence of excitatory events generated by a homogeneous stochastic process. We could show that the ISI distributions are not compatible with the assumption that the intervals between excitatory events are exponentially distributed, i.e., that the excitatory events are generated by a homogeneous (or a doubly stochastic) Poisson process, at least when the refractory period is modeled as consisting of an absolute refractory period followed by a relative refractory period during which the probability of a spike occurring given an excitatory event increases exponentially from 0 to 1. Such a refractory function has been implemented in numerous models and is consistent with experimental data (see references in results). Theoretically, however, our finding that the reduced version of the gamma-mixture model provides a much better fit than the full version of the Poisson model could be compatible with the assumption that the excitatory events are generated by a homogeneous (or a doubly stochastic) Poisson process. In this case, a relative refractoriness is required that never removes more than two events in a row from the Poisson process. Such a refractory behavior has, to the best of our knowledge, never been observed in real neurons. We therefore favor the alternative interpretation of our findings, viz. that the excitatory events are generated by a homogeneous stochastic process where the distribution of their intervals is a mixture of an exponential and a gamma distribution with shape factor 2, both with the same scaling parameter λ_{E} and that these excitatory events then interact with a conventional refractory period.

As explained before, this gamma-mixture of interval distributions can emerge from a primary homogeneous Poisson point process with rate λ_{E}, if some of these primary events are annihilated by a subsequent process or fail to trigger secondary events further down a cascade of precursor steps that might be required to produce the ultimate excitatory events. However, primary events cannot be annihilated at random because the resulting distribution of secondary events would still be exponential, merely with a smaller scaling factor. Instead, primary events must be removed obeying the “never-two-in-a-row” rule in order for the distribution of the intervals between the remaining events to be a mixture of a gamma distribution with β = 2 and an exponential distribution of common scaling factor. Apart from effects produced by a recording artifact and by preferred intervals, the ISI distributions were well, although not fully, explained by, this model. This gamma-mixture model was very successful for describing the ISI distributions of cat AN fibers (Heil et al. 2007). Inspection of ISI histograms or probability density functions obtained from the spontaneous activity of AN fibers in the guinea pig (Manley and Robertson 1976; Figs. 2, 5, and 6), as well as the pigeon (Gummer 1991b; Fig. 2), suggests that the gamma-mixture model would also be successful in other mammalian and avian species.

### Possible physiological basis of the excitatory process of the gamma-mixture model

It should be clear that our model is not a complete description of the processes of spike generation in primary afferents. We only attempted to model a single aspect of their spiking behavior, viz. the distribution of ISIs during spontaneous activity. Nevertheless, it is interesting to speculate about some of the underlying physiological mechanisms. In the gamma-mixture model, each ultimate excitatory event triggers a spike in the AN fiber unless the fiber is refractory. Therefore such excitatory events can be conceived of as excitatory postsynaptic potentials (EPSPs) of amplitudes large enough to trigger a spike in the AN fiber when it is not refractory. For simplicity, we will refer to these EPSPs as suprathreshold (even if they fail to trigger a spike during the refractory period). The question thus arises: how could such a mixture of distributions be generated in the system providing the suprathreshold EPSPs in AN fibers? In principle, the annihilation of some of the primary events, in the sense above and according to the never-two-in-a-row rule, could occur both postsynaptically or presynaptically. In the first scenario, the primary events would be transmitter release events from presynaptic sites in the hair cell and their intervals would be exponentially distributed with scaling factor λ_{E}. It should be stressed that a transmitter release event in this sense is not necessarily identical to the release of a single quantum of transmitter. Release events could be uniquantal but could also be instances where the transmitter content of two or more vesicles is released simultaneously. In mouse inner hair cells, in bullfrog saccular hair cells as well as in rat rod bipolar cells, there is evidence for such multivesicular, coincident, or coordinated release (Glowatzki and Fuchs 2002; Keen and Hudspeth 2006; Neef et al. 2007; Singer et al. 2004). In this scenario, a fraction *b* of these primary (transmitter release) events would have to be annihilated such that only the remaining secondary events (suprathreshold EPSPs) will lead to spikes in the afferent fiber, unless they fall into the refractory period that follows each spike. Several possibilities exist for how this could be achieved, but in light of the fact that the annihilation mechanism has to obey the never-two-in-a-row rule, they all appear unlikely. *1*) Annihilation could be achieved by rapid removal of transmitter from the synaptic cleft, but it is difficult to conceive of a removal mechanism that could obey the rule. *2*) Each transmitter release event might lead to an EPSP, sub- or suprathreshold, so that the EPSP intervals were also exponentially distributed (as has been suggested in some preparations; Fatt and Katz 1952; Holt et al. 2006). Now, it would have to be assumed that a fraction *b* of these EPSPs is too small to trigger a spike even if the AN fiber is not refractory. Because EPSPs likely differ in amplitude, because of variation in the amount of transmitter released (e.g., because of variation in the number of simultaneously released vesicles) or variation in its effectiveness (e.g., because of accumulation of transmitter or receptor desensitization; Holt et al. 2006), some of them may indeed be too small to trigger a spike. However, the never-two-in-a-row rule requires that subthreshold EPSPs never occur in succession, and it is very difficult to see how this could be ensured. *3*) It is also conceivable that the postsynaptic spike generator might be responsible for the gamma-mixture. A gamma distribution with β = 2 can emerge from an exponential distribution via a perfect integrator that needs to integrate over *n* = 2 input events before reaching threshold and being reset to zero (Koch 1999; Stein 1965). Threshold fluctuations between *n* ≤ 1 and *n* ≤ 2 of such an integrator would produce a mixture of gamma (β = 2) and exponentially distributed intervals between its output events with a common scaling factor. However, this scenario requires the input events to have negligible decay, whereas excitatory postsynaptic currents and potentials in the peripheral dendrites of primary afferents decay very rapidly (Glowatzki and Fuchs 2002; Goutman and Glowatzki 2007; Holt et al. 2006; Keen and Hudspeth 2006). Therefore fluctuations in the threshold for spike initiation (see e.g., Bruce et al. 1999b) are also unlikely to be the physiological mechanism.

These considerations favor the second scenario, namely that the distribution of intervals between transmitter release events is already a mixture of a gamma and an exponential distribution and so is the distribution of intervals between the resulting EPSPs. Transmitter-containing vesicles undergo a series of steps, including trafficking to the active zone, docking with the cell's membrane, and a priming step, before they finally fuse with the cell membrane and release their transmitter (Prescott and Zenisek 2005). The intervals between events on an early (primary) level in such a cascade might have Poisson statistics. For example, Prescott and Zenisek (2005) have suggested that, at ribbon synapses, the replenishment of vesicles at the ribbon occurs through Brownian vesicle motion leading to random collision of vesicles with the ribbon. If the transitions from such a primary level to the next in the maturational cascade of transmitter release were to fail with probability *b*, but happen the next time an event at the primary level occurred, the distribution of the intervals between the secondary events would automatically be a mixture of a gamma distribution with β = 2 and fraction *b* and an exponential distribution with fraction 1 − *b*, where both have the same scaling factor. We believe it is unlikely that the gamma component is caused by activity of efferent neurons projecting with thin axons from the brain stem to the basilar papilla and terminating on tall hair cells, although essentially nothing is known about their physiology (Code 1997; Raabe and Köppl 2003).

### Proportion of gamma-distributed intervals

In the cat, where the majority of AN fibers is driven by a single ribbon synapse each (Liberman et al. 1990), we provided some evidence that the fraction *b* of gamma-distributed intervals between excitatory events is constant across AN fibers, i.e., independent of spontaneous discharge rate or CF (Heil et al. 2007). This suggested to us that all ribbon synapses in cat inner hair cells might operate in a similar manner, just at different rates. The mean *b*, which we obtained from our sample of cat AN fibers, was 0.40 ± 0.17 and the median was 0.43 (Heil et al. 2007). The fraction *b* of gamma-distributed intervals obtained here from barn owl AN fibers was considerably smaller (viz., 0.202 ± 0.140 from fits to individual CDFs without bias correction, ∼0.225 with bias correction, 0.224 from the mean CDF, and 0.238 from the CV-skew analysis). Furthermore, in the barn owl, the estimates of *b* were positively correlated with the mean ISI, reflecting a tendency for *b* to be smaller the higher the spontaneous rate, and with CF, reflecting the general trend toward lower spontaneous rates with increasing CF. A simple explanation for these differences might be that the ribbon synapses in the two species function differently, with additional differences as a function of spontaneous rate and CF in the barn owl. However, there is no experimental evidence to suggest that ribbon synapses in hair cells of different sources differ in fundamental properties (reviewed by Glowatzki et al. 2008). A more parsimonious assumption is therefore that they operate in a similar manner, i.e., with a similarly large *b*. If ribbon synapses in cat and barn owl operated in a similar manner, say with the value of *b* = 0.43 found to be appropriate for capturing the ISI distributions of AN fibers in the cat (Heil et al. 2007), other mechanisms are required to account for the smaller estimates of *b* and for their correlation with mean ISI and CF in the barn owl. Several mechanisms, which are not mutually exclusive, are conceivable.

A reduction from a hypothetical, large value of *b* at the level of transmitter release events from a single ribbon synapse to its estimate derived from ISIs can result from random annihilation of events at any level in between, for example, by the processes discussed above (random transmitter removal, random fluctuations of EPSP amplitudes, spike threshold fluctuations, or efferent activity). The observed correlation of *b* estimates with mean ISI in the barn owl would require that the proportion of random annihilations systematically increases with increasing spontaneous rate. From an energetic point of view, however, annihilations at these levels seem suboptimal. A second possibility to account for differences in *b* estimates between cat and owl, under the assumption that individual ribbon synapses operate in a similar manner, arises if single AN fibers in the barn owl were driven by more than one synaptic ribbon. If two or more independent ribbons, each operating with a *b* of say 0.4–0.5, would drive a single primary afferent, the estimate of *b* that is derived from the ISI distribution will be lower than the value of *b* of the individual ribbons. The more ribbons, the lower *b*. In other words, the ISI distributions will be more compatible with an exponential distribution, generated by a homogeneous Poisson process, and modified by refractoriness. In the chicken, this Poisson model seems to provide a satisfactory description of AN fiber activity (Avissar et al. 2007).

Quantitative data on afferent innervation density and synaptic ribbon numbers in the chicken (*Gallus gallus*) suggest that this could indeed be the case. Here, each tall hair cell on the inner (neural) side of the basilar papilla receives only one to three afferent terminals (Fischer 1992), but contains ∼15 presynaptic dense bodies (Martinez-Dunst et al. 1997). This implies 5–15 ribbons per primary afferent. Less quantitative observations on other bird species support the notion that this may be typical for birds. Several papers explicitly mention, or show micrographs that illustrate, that more than one ribbon can be associated with the same afferent ending (Chandler 1984; Fischer 1992; Fischer et al. 1992; Takasaka and Smith 1971). In addition, a minority of fibers typically contacts more than one hair cell in the apical region of the basilar papilla (Fischer 1992, 1994, 1998; Fischer et al. 1992). Each of those fibers would likely be excited by more than one synapse.

In light of the above considerations, the correlation of *b* with mean ISI (Fig. 6*B*) and with CF could indicate an increase in the number of ribbons per primary afferent from the base of the basilar papilla to its apex. Such a gradient could arise via a decrease in the number of afferents per tall hair cell along this axis, even if the number of ribbons per tall hair cell were constant. In this context, it is noteworthy that perhaps the most striking and unusual feature of afferent innervation in the barn owl is the higher number of afferent terminals on tall hair cells, compared with other bird species, especially in regions of the papilla corresponding to CFs of ∼5–8 kHz. Here, ≤20 terminals were observed on the most neural hair cells (Fischer 1994).

Differences in the number of ribbons per primary afferent might also contribute to differences in spontaneous rates in birds, but other factors likely dominate (Heil and Neubauer 2001; Meddis 2006; Saunders et al. 2002 for a discussion), as highlighted by mammals, where only a single ribbon drives most fibers (Liberman et al. 1990), yet there are large differences in spontaneous rates among primary afferents (Joris and Yin 1992; Kiang et al. 1965; Liberman 1978; Relkin and Doucet 1991; Taberner and Liberman 2005; Winter and Palmer 1991; Yates 1991).

### Refractoriness

The shortest ISIs observed in these data averaged to 0.68 ± 0.11 ms. In starling AN fibers, the same measure was reported to vary between 0.5 and 2.2 ms, with the bulk of the measures falling between 0.8 and 1.6 ms (Manley et al. 1985). An average was not provided, but would be >1 ms, as judged from Figs. 8*B* and 9*B* in that study. In the pigeon, the shortest ISIs averaged 0.80 ± 0.06 and 0.97 ± 0.06 ms for fibers without and with preferred intervals, respectively (Gummer 1991a). In all those cases, the shortest ISIs can only constitute a conservative estimate of the upper limit of the absolute refractory period, because of the recording artifact. This artifact also prevented us from fitting the models to very short ISIs. Therefore we were not able to estimate the durations of the absolute and relative refractory periods individually. However, we could estimate the sum of the absolute refractory period and of the mean duration of the relative refractory period with some confidence. Across the 136 AN fibers we obtained a mean estimate of 0.365 ± 0.266 ms before bias correction, corresponding to ∼0.4 ms after bias correction, an estimate of 0.535 ms from the mean CDF, and 0.543 ms from the CV-skew analysis. The latter two estimates are probably most reliable. They are remarkably brief. In the cat, for example, we estimated a mean refractory period of ∼1.24 ms, with an absolute refractory period of 0.59 ms and a mean duration of the relative refractory period of 0.65 ms (Heil et al. 2007), in line with estimates from other studies of mammalian species, including humans (Brown 1994; Cartee et al. 2000; Dynes 1996; Miller et al. 2001; Morsnowski et al. 2006; Parkins 1989; Shepherd et al. 2004). Similar numbers have also been used in some modeling studies (Bruce et al. 1999a; Litvak et al. 2003; Lütkenhöner et al. 1980; Meddis 2006; Meddis and O'Mard 2005; Schroeder and Hall 1974). Hence, the estimated refractory period in the barn owl is two to three times shorter than in the cat. This is remarkable, particularly because spike durations of barn owl and cat AN fibers seem to be similar. From a sample of 17 barn owl AN fibers from a population recorded with comparable methods to this study (Köppl, unpublished data), we derived a median width at half height of 190 μs. From Fig. 14 of Gaumond et al. (1982), we derived a corresponding estimate of 190 μs for the cat AN fiber action potential illustrated there. The physiological mechanisms enabling the brief refractory period in the barn owl are unknown, but given the similar durations of action potentials in barn owl and cat, they are unlikely to be caused by mechanisms that might also affect the duration of the action potential. In any event, the brevity of the refractory period in the barn owl may reflect an adaptation to the needs for fast temporal processing, as evident in its phase locking to high frequencies (Köppl 1997c).

### Recording artifact affects short ISIs

We observed unnaturally steep, often step-like, increases of the CDF at or around the shortest ISIs observed in a given sample of spontaneous activity (Fig. 2*A*). These rapid increases manifest themselves as an “early peak” (Köppl 1997a) in the ISI histogram and as a sharp nonmonotonic peak in the hazard function (e.g., Fig 2*B*). We argued that these effects are artifacts related to the recording technique. Independent evidence for the presence of such a recording artifact was provided by analyses of the oscillations of the hazard rates of some AN fibers (Fig. 2, *C–F*), because they are most likely caused by spontaneous membrane-potential oscillations associated with an electrical tuning mechanism in the hair cells (Köppl 1997a). When a fixed trigger threshold is used, as has been and still is done in numerous studies, including Köppl (1997a), the second of a pair of consecutive spikes is triggered at a later phase than the first one, if the second spike has a smaller amplitude than the first one and/or if it “rides” on the afterpotential following the first spike (see also Gaumond et al. 1982). Therefore the recorded ISI will be longer than the true ISI. Because spike amplitude is decreased at short ISIs (Gaumond et al. 1982; Shepherd et al. 2004; Siegel 1992) and the afterpotential is of finite duration, this effect will be most prominent at very short ISIs. The use of micropipettes with rather high impedances (of 50–100 MΩ by Köppl 1997a) might intensify the artifact, because their long time constants could prolong the apparent afterpotential. In support of this conclusion, it is noteworthy that studies in mammals and birds, in which an initial sharp nonmonotonic peak in the hazard functions or in ISI histograms has also been reported (Gaumond et al. 1982; Gray 1967; Klinke et al. 1994; Li and Young 1993; Miller and Wang 1993; Prijs et al. 1993; Young and Barta 1986), have all used high-impedance micropipettes and a fixed trigger threshold to record the fiber activity and the timing of spikes. In our previous study of cat AN fibers (Heil et al. 2007), we also used a fixed trigger threshold but used tungsten electrodes of relatively low impedance (4–7 MΩ). In that study, the artifact was relatively small compared with the one seen here in the barn owl or to the ones published in the literature cited above. A second factor that may contribute to the apparent delay of the second spike at short ISIs are possible differences in conduction velocity, with that for the second spike being lower, as discussed by Gaumond et al. 1982. However, the presence of only relatively small apparent delays in our cat study (Heil et al. 2007) argues against possible differences in conduction velocity as a major contributor. In any event, it is also conceivable that the artifact and the use of high-impedance pipettes with long time constants contribute to the differences between experimental and simulated data in residuals (Fig. 8) and in the course of the SD of the logarithms of (ISI–*t*_{RP}) as functions of the ISI (Fig. 11). The best solution to avoid the artifact in future studies would be to record and store the entire electrode signal and to define the ISIs as the intervals between the peaks of the spike waveforms.

### Conclusions

We showed here that the ISI distributions from spontaneous activity of AN fibers of the barn owl, even of those not showing preferred intervals, are not compatible with the assumption that excitatory events are generated by a homogeneous Poisson point process and that lead to spikes unless the fiber is refractory. The ISI distributions are better, although not fully, explained as resulting from the interaction of a brief refractory period (∼0.5 ms) and excitatory events generated by a homogeneous stochastic process where the distribution of interevent intervals is a mixture of an exponential and a gamma distribution with shape factor 2, both with the same scaling parameter. An artifact, caused by the common technique of recording with high-impedance micropipettes in combination with a Schmitt trigger, severely affected short ISIs and precluded a more detailed analysis. For cat AN fibers, the same model provided an accurate description of the ISI distributions (Heil et al. 2007). Furthermore, the mean proportions of exponentially versus gamma distributed intervals in the mixture were different for cat and owl, and they were constant across fibers in the cat while they covaried with mean spontaneous rate and with characteristic frequency in the owl. Future studies may shed light onto the origin of these differences, which could reside, for example, in differences in the number of ribbon synapses providing excitation to cat versus owl primary afferents, as well as onto the molecular or cellular mechanisms that give rise to the gamma-mixture.

## Appendix A

In this appendix we show that, provided the mean relative refractory period is short compared with the mean interval between excitatory events, i.e., 1/λ_{R} ≪ 1/λ_{E}, or equivalently λ_{R} ≫ λ_{E}, *Eq. 3*, for ISIs with durations *t* ≫ 1/λ_{R}, can be well approximated by an exponential distribution with a deadtime of *t*_{RP} ≅ *t*_{D} + 1/λ_{R}.

The general-gamma distribution, described by *Eq. 3*, is reproduced here for convenience (A1) The part of the equation that is valid for *t* ≥ *t*_{D} can be rewritten as (A2) As (*t* − *t*_{D}) increases, the middle summand in the right side of *Eq. A2* approaches zero much faster than the last summand if λ_{R} ≫ λ_{E}, so the equation simplifies to (A3) The factor λ_{R}/(λ_{R} − λ_{E}) is equivalent to an extension of the deadtime as shown next. *Equation A3* can be rewritten as (A4) or as (A5) The term ln[λ_{R}/(λ_{R} − λ_{E})] in the exponent of *Eq. A5* can be rewritten as ln[1 + λ_{E}/(λ_{R} − λ_{E})]. When λ_{R} ≫ λ_{E}, then λ_{E}/(λ_{R} − λ_{E}) ≪ 1, and ln[1 + λ_{E}/(λ_{R} − λ_{E})] can be well approximated by the first term of a Taylor series, viz. λ_{E}/(λ_{R} − λ_{E}), which can be further approximated by λ_{E}/λ_{R}. With these approximations, *Eq. A5* simplifies to (A6) Hence, *Eq. A1* can be approximated by (A7) Thus a fit of *Eq. A7* restricted to the upper part of an ISI distribution allows the retrieval of the approximate mean refractory period, because *t*_{RP} ≅ *t*_{D} + 1/λ_{R}, but not of the individual durations of its two components.

### Appendix B

In this appendix we show that, provided the mean relative refractory period is short compared with the mean interval between excitatory events, i.e., 1/λ_{R} ≪ 1/λ_{E}, or equivalently λ_{R} ≫ λ_{E}, *Eq. 7*, for ISIs with durations *t* ≫ 1/λ_{R}, can be well approximated by *Eq. 8* with a deadtime of *t*_{RP} ≅ *t*_{D} + 1/λ_{R}.

*Equation 7* is reproduced here for convenience (B1) The first summand, with the weighting factor of (1 − *b*), describes the distribution of the fraction of ISIs that results from exponentially distributed intervals between excitatory events (cf. *Eq. 3*). For ISIs ≫ 1/λ_{R} and λ_{R} ≫ λ_{E}, it simplifies to the exponential distribution, given by *Eq. A7*, as detailed in appendix A. The second summand, with the weighting factor of *b*, describes the distribution of the fraction of ISIs that results from the gamma distributed intervals between excitatory events. It is composed of three subsummands, viz., S1 = 1, S2, and S3. The middle subsummand, S2, can be written as It is easy to see, again by comparison with appendix A, that for (*t* − *t*_{D}) ≫ 1/λ_{R} and λ_{R} ≫ λ_{E}, S2 can be approximated by (B2) The third subsummand, S3, can be rewritten as and approximated, following the procedures explained in appendix A, by (B3) The sum of the approximations of S2 and S3 is thus given by which, for λ_{R} ≫ λ_{E}, further approximates to (B4) Substituting the corresponding terms in *Eq. B1* by these approximations yields (B5) Thus a fit of *Eq. B5* restricted to the upper part of an ISI distribution allows the retrieval of the approximate mean refractory period, because *t*_{RP} ≅ *t*_{D} + 1/λ_{R}, but not of the individual durations of its two components. In addition, the fraction *b* of gamma-distributed intervals can also be derived correctly from that portion of the ISI distribution.

## GRANTS

This work was supported by Deutsche Forschungsgemeinschaft Grant SFB-TR31 A6 to P. Heil. C. Köppl conducted the experiments and provided the raw data, H. Neubauer developed the model, the curve fitting procedure, and many of the analysis tools. H. Neubauer and P. Heil analyzed the data. P.Heil, H. Neubauer, and C. Köppl wrote the paper.

## Acknowledgments

We thank Prof. Dr. Georg Klump for help in recovering the spike times and for excellent custom hardware and software used during the experiments and two anonymous reviewers for thorough reviews of, and constructive comments on, two earlier versions of the manuscript.

- Copyright © 2009 the American Physiological Society

## REFERENCES

- Avissar et al. 2007.↵
- Ayyub 1997.↵
- Bartlett 1949.↵
- Bi 1989.↵
- Brown 1994.↵
- Bruce et al. 1999a.↵
- Bruce et al. 1999b.↵
- Carney 1993.↵
- Cartee et al. 2000.↵
- Chandler 1984.↵
- Code 1997.↵
- Cox 1962.↵
- Delgutte 1996.↵
- Dynes 1996.↵
- Fatt and Katz 1952.↵
- Fischer 1992.↵
- Fischer 1994.↵
- Fischer 1998.↵
- Fischer et al. 1992.↵
- Fuchs 2005.↵
- Gaumond et al. 1983.↵
- Gaumond et al. 1982.↵
- Geisler 1981.↵
- Geisler 1998.↵
- Geisler et al. 1985.↵
- Glowatzki and Fuchs 2002.↵
- Glowatzki et al. 2008.↵
- Goutman and Glowatzki 2007.↵
- Gray 1967.↵
- Gummer 1991a.↵
- Gummer 1991b.↵
- Heil and Neubauer 2001.↵
- Heil et al. 2007.↵
- Holt et al. 2006.↵
- Jackson and Carney 2005.↵
- Javel 1986.↵
- Joris and Yin 1992.↵
- Keen and Hudspeth 2006.↵
- Kiang et al. 1965.↵
- Klinke et al. 1994.↵
- Koch 1999.↵
- Köppl 1997a.↵
- Köppl 1997b.↵
- Köppl 1997c.↵
- Köppl et al. 1993.↵
- Köppl and Yates 1999.
- Krishna 2002.↵
- Kuhlmann et al. 2002.↵
- Lehmann 2002.↵
- Li and Young 1993.↵
- Liberman 1978.↵
- Liberman et al. 1990.↵
- Liberman and Kiang 1978.↵
- Litvak et al. 2003.↵
- Lowen and Teich 1991.↵
- Lowen and Teich 1992.↵
- Lütkenhöner et al. 1980.↵
- Manley et al. 1985.↵
- Manley and Köppl 1998.↵
- Manley and Robertson 1976.↵
- Martinez-Dunst et al. 1997.↵
- Meddis 2006.↵
- Meddis and O'Mard 2005.↵
- Miller et al. 2001.↵
- Miller and Wang 1993.↵
- Molnar and Pfeiffer 1968.↵
- Morsnowski et al. 2006.↵
- Moser et al. 2006.↵
- Moser et al. 2006.
- Mountain and Hubbard 1996.↵
- Neef et al. 2007.↵
- Nouvian et al. 2006.↵
- Parkins 1989.↵
- Perkel et al. 1967.↵
- Prescott and Zenisek 2005.↵
- Prijs et al. 1993.↵
- Raabe and Köppl 2003.↵
- Relkin and Doucet 1991.↵
- Richter et al. 1995.↵
- Sachs 2002.↵
- Saunders et al. 2002.↵
- Schmich and Miller 1997.↵
- Schoonhoven et al. 1997.↵
- Schroeder and Hall 1974.↵
- Shepherd et al. 2004.↵
- Siegel 1992.↵
- Singer et al. 2004.↵
- Stein 1965.↵
- Stephens 1986.↵
- Sterling and Matthews 2005.↵
- Taberner and Liberman 2005.↵
- Takasaka and Smith 1971.↵
- Tanaka and Smith 1978.↵
- Teich et al. 1990.↵
- Winter and Palmer 1991.↵
- Yates 1991.↵
- Young and Barta 1986.↵
- Zeddies and Siegel 2004.
- Zhang et al. 2001.↵