## Abstract

Reliable accounts of the variability observed in neural spike trains are a prerequisite for the proper interpretation of neural dynamics and coding principles. Models that accurately describe neural variability over a wide range of stimulation and response patterns are therefore highly desirable, especially if they can explain this variability in terms of basic neural observables and parameters such as firing rate and refractory period. In this work, we analyze the response variability recorded in vivo from locust auditory receptor neurons under acoustic stimulation. In agreement with results from other systems, our data suggest that neural refractoriness has a strong influence on spike-train variability. We therefore explore a stochastic model of spike generation that includes refractoriness through a recovery function. Because our experimental data are consistent with a renewal process, the recovery function can be derived from a single interspike-interval histogram obtained under constant stimulation. The resulting description yields quantitatively accurate predictions of the response variability over the whole range of firing rates for constant-intensity as well as amplitude-modulated sound stimuli. Model parameters obtained from constant stimulation can be used to predict the variability in response to dynamic stimuli. These results demonstrate that key ingredients of the stochastic response dynamics of a sensory neuron are faithfully captured by a simple stochastic model framework.

## INTRODUCTION

Sensory neurons provide an animal with the primary representation of its environment and internal state. Any computation and behavioral decision has to be based on this representation. Yet even under controlled laboratory conditions, repeated stimulus presentations often lead to a considerable trial-to-trial variability of neural responses. This is particularly true for cortical neurons, which typically discharge with high variability under in vivo conditions (Buracas et al. 1998; Holt et al. 1996; Shadlen and Newsome 1998; Stevens and Zador 1998). As expected for Poisson-like random processes, coefficients of variation near unity and spike-count variances equal to or even above the mean count have been found (Softky and Koch 1993; Teich et al. 1996).

On the other hand, neural response characteristics depend sensitively on the temporal features of the stimulus pattern (Mainen and Sejnowski 1995), which is of special importance for sensory neurons receiving highly structured dynamic stimuli (de Ruyter van Steveninck et al. 1997; Kara et al. 2000; Machens et al. 2001; Warzecha et al. 2000). Indeed, for appropriate inputs, individual spikes can be highly reliable and precisely timed (Berry et al. 1997; Reinagel and Reid 2002), resulting in spike-count variances far below the mean spike count (de Ruyter van Steveninck et al. 1997; Warzecha and Egelhaaf 1999).

A central question in understanding neural coding principles therefore concerns the nature, origin, and computational implications of neural variability. To investigate these aspects, it is essential to advance stochastic descriptions of neural response dynamics that are valid over a broad spectrum of conditions including different neural activation strengths and different types of temporal stimulus modulation.

Based on in vivo recordings from locust auditory receptor neurons, a model system for the auditory periphery of insects (Michelsen 1971; Römer 1976; Ronacher and Krahe 2000; Stumpner and von Helversen 2001), we systematically explore a phenomenological description that captures the spike-train statistics at different average firing rates and that applies to constant as well as temporally modulated sound stimuli. Because the receptors are the first stage of the auditory system, the experimentally observed response variability is not inherited from other neurons, but must be caused by intrinsic processes. Furthermore, the accuracy and reliability of responses in the investigated cells has been shown to be strongly influenced by the stimulus statistics (Machens et al. 2001, 2003), which poses a particular challenge for the mathematical description.

To characterize the receptor dynamics within a general theoretical framework, we first assess the spike-train variability in response to constant-intensity stimuli over a large range of sound frequencies and intensities. The observed interspike interval (ISI) statistics suggests that spike generation can be modeled by a renewal process (Cox 1962). To account for neural refractoriness, recovery functions (Berry and Meister 1998; Johnson 1996) are incorporated into the framework. We use a particularly simple realization where for each cell, one recovery function is determined from a single ISI distribution that is obtained from stimulation with a constant sine tone. As shown by our data, a renewal process based on such a recovery function faithfully describes the shape of ISI distributions for arbitrary sound frequencies and intensities. Combining this stochastic spike generator with a deterministic stimulus encoder allows us to calibrate the model neurons with independent measurements of the receptors' input–output relation. Our results are therefore direct predictions from the stimulus and do not rely on an observed poststimulus time histogram (PSTH). The general model accurately accounts for spike-train variability in response to a variety of both constant and dynamic stimuli.

## METHODS

### Electrophysiology

All experiments were performed on adult locusts (*L. migratoria*). Legs, wings, head, and gut were removed to immobilize the animals and to facilitate access to the metathoracic ganglion and auditory nerve. Preparations were fixed with wax, ventral side down, onto a Peltier element.

The dorsal part of the thorax was opened to expose the metathoracic ganglion and the auditory nerve, which was fixed with a custom-made forceps mounted on a micromanipulator. During the experiments, preparations were kept at a fixed temperature of about 30°C. Acoustic stimuli were presented by loudspeakers [Esotec D-260, Dynaudio (Skanderborg, Denmark) on a DCA 450 amplifier (Denon Electronic GmbH, Ratingen, Germany)]. Receptor cells were recorded intracellularly from the axon in the auditory nerve with standard glass microelectrodes (borosilicate, GC100F10; Harvard Apparatus, Edenbridge, UK), filled with a 1 M KCl solution (30–60 MΩ resistance). Neural responses were amplified (BRAMP-01, NPI Electronis, Tamm, Germany) and recorded by a data-acquisition board (National Instruments, PCI-MIO-16E-1) with a sampling rate of 10 kHz. Stimulus generation, spike detection, and data analysis were performed using custom-made software. Spike times were determined with a temporal resolution of 0.1 ms and stored for off-line analysis. The response latency, estimated separately for each neuron from the shortest observed time difference between stimulus onset and the first spike, was subtracted from the spike times. In all experiments, the first 300 ms of the response were discarded to minimize the influence of firing-rate adaptation (Benda 2002). The experimental protocol complied with German law governing animal care.

### Stimulus design

Pure tones whose frequencies ranged from 1 to 40 kHz were used to determine threshold curves and thus to classify each cell as one of 4 standard receptor cell types (Römer 1976). Types I–III have their greatest sensitivity in the lower-frequency range between 3 and 8 kHz with differences in absolute sensitivity and the exact location of the sensitivity maximum; type IV is most sensitive at high frequencies around 15 kHz.

In the first set of experiments, stimuli were pure tones of constant intensity and 1-s duration, followed by a 1-s–long quiet pause. One group of 8 cells was stimulated with at least 4 different sound frequencies (between 3 and 11 kHz for low-frequency receptors; between 12 and 25 kHz for high-frequency receptors). These stimuli were repeated up to 10 times. For another group of 11 cells, the sound frequency was set to the characteristic frequency (i.e., the frequency of highest sensitivity). In this paradigm, stimuli were repeated up to 50 times. For both groups, 5 to 10 different sound intensities were used to cover the whole range of firing rates from threshold to saturation. From the recorded spike trains, ISI distributions and coefficients of variation (see *Data analysis*) were calculated.

In a second set of experiments, stimuli were pure tones at the cell's characteristic frequency whose amplitudes were modulated by Gaussian white noise with a cutoff frequency of 400 Hz. The standard deviation of the noise signal around the mean intensity *I*_{0} was either 3 dB (5 cells) or 5 dB (one cell). To restrict the amplitudes within a finite range, the tails of the Gaussian distributions were cut off at 4 SDs. The modulation depths (defined as the central 95% of the amplitude distribution) of the stimuli were 11.76 dB [3 dB SD stimulus] and 19.6 dB (5 dB SD stimulus). Stimuli were again 1 s long, separated by 1-s–long pauses, and now repeated 60 times.

### Data analysis

##### FIRING RATE AND SPIKE-TRAIN VARIABILITY.

Firing rates were calculated as trial averages of spike counts in sliding windows of 10 ms length. The variability of interspike intervals (ISIs) was analyzed using the coefficient of variation (*CV*), which is defined as the ratio between the standard deviation σ_{ISI} of the ISI distribution and its average 〈*ISI*〉 (1) To quantify the spike-count variability, Fano factors (Fano 1947) were calculated. The Fano factor *F*(*T*) is derived from the spike-count distribution for a certain counting time *T* by dividing the variance of the spike count σ_{N(T)}^{2} by its average 〈*N*(*T*)〉 (2) The errors of the measured *CV* values and Fano factors were determined using the bootstrap method. New synthetic data sets, numbering 100, were constructed from the original data (ISIs or spike counts) by drawing randomly with replacement. Each synthetic data set had the same size as the original data. *CVs* or Fano factors were calculated for each synthetic data set, and the SDs of these reevaluated quantities were taken as an estimate of the error of the original *CV*/Fano factor.

To assess the quality of the model's Fano factor predictions, we quantified the absolute deviation of the model Fano factor *F*_{m}(*t*) from the experimental Fano factor *F*_{e}(*t*) for each time point *t*, and scaled this deviation by the error of the experimental Fano factor *Error* [*F*_{e}(*t*)] at the same point in time. The scaled deviation was then averaged over time to yield the mean relative prediction error (3) To quantify spike timing reliability in experimental and model spike trains, we used a correlation-based measure introduced by Schreiber et al. (2003): The spike trains obtained from *N* repeated presentations of a stimulus are convolved with a Gaussian filter with width (SD of the Gaussian) σ_{c}, yielding a smoothed representation *s*_{i}(*t*) of the spike trains, where the index i enumerates the repeated stimulus representations. These representations can also be viewed as vectors *s*i, where every component corresponds to a point in time. Subsequently, the inner product is taken between all possible pairs of spike trains *s*i, *s*j and each inner product is then divided by the norms of the 2 spike trains of the representative pair. The average over all pairs of spike trains is then taken as the reliability *R*_{corr} (4) This measure takes values between 0 (minimum reliability) and 1 (maximum reliability). We used values of 0.5 to 3 ms for the filter width σ_{c} to quantify the spike-timing reliability on different timescales.

##### INDEPENDENCE OF INTERSPIKE INTERVALS.

Modeling stochastic neural responses is greatly facilitated if the ISIs can be assumed to be independent under constant stimulation (see *Model*). Two tests are particularly suited to investigate whether this assumption is satisfied by the spike trains under study. In the first test, correlations between an interspike interval *ISI*_{k} and any of its successors *ISI*_{k+j} are analyzed by calculating serial correlation coefficients *r*_{j} for lags *j* up to 20 (*j* ≥ 0) (5) where σ_{ISI}^{2} denotes the variance of the ISIs, and 〈 · 〉 stands for the average over *k*. Because independence implies the lack of correlations, serial correlation coefficients for *j* ≥ 1 that significantly differ from zero indicate that the ISIs are not independently distributed.

For a renewal process with independent ISIs, the autocorrelation function can be calculated from the ISI distribution (Perkel et al. 1967). Therefore we used a comparison of the autocorrelation function *C*(Δ) of the observed spike train (mean not subtracted) to the autocorrelation function *C*_{iid}(Δ) of a spike train with the same ISI distribution and independent, identically distributed (iid) ISIs as a second test for independence of ISIs. *C*_{iid}(Δ) is equal to the iterative convolution of the ISI distribution *P*_{ISI}(Δ) (see Perkel et al. 1967 for a detailed derivation) (6) In practice, a finite number of terms is sufficient because neurons do not exhibit infinitely small ISIs. If *C*(Δ) and *C*_{iid}(Δ) differ significantly, the observed ISIs cannot be assumed independent.

### Model

Our model is based on the assumption that spikes are generated stochastically at time *t* with a probability that depends on the product of 2 terms only: the strength *q*(*t*) of an “effective stimulus” at that very moment and a term that depends on the length of the interval Δ since the last spike (generated at time *t*_{last} = *t* − Δ). Earlier spikes or the stimulus strength between *t*_{last} and *t* do not influence spike generation at time *t*. For time-independent stimuli [i.e., *q*(*t*) = *q*], spike generation is thus a “renewal process” (Cox 1962). This implies that ISIs are independent. For time-dependent stimuli, the latter property is no longer true and one rather speaks of a “modulated renewal process” (Reich et al. 1998).

The functional dependency of the spike probability on Δ is denoted by *w*(Δ). This memory term, or “recovery function” (Berry and Meister 1998; for a comparison with the “hazard function,” see Gerstner and Kistler 2002; Johnson 1996), captures the influence of refractoriness on the generation of the next action potential. Mathematically, spike generation is thus described by a probability per unit time (the “hazard”) ρ(*t*|*t*_{last}) that is conditional on the last spike occurring at time *t*_{last} (7) The values of the recovery function *w* range between zero and unity. Within the absolute refractory period τ_{a}, even arbitrarily large stimuli cannot elicit a spike. To capture this property, *w*(Δ) vanishes for 0 ≤ Δ ≤ τ_{a}. It then rises monotonically to account for the relative refractory period during which the neuron relaxes back to its normal level of excitability. Note that the probability of spike generation explicitly depends on the time that has passed since the last spike. Unlike for a Poisson process, individual spikes are therefore not independent. However, ISIs still are independent under constant stimulation because refractoriness is not cumulative over multiple spikes.

For a renewal process, it is possible to compute the recovery function directly from the ISI distribution *P*_{ISI}(Δ) obtained under constant stimulation [*q*(*t*) = *q*], as has been discussed in the literature (Berry and Meister 1998; Gerstner and Kistler 2002; Johnson 1996) (8) *Equation 8* can be inverted to calculate the ISI distribution from the recovery function (9) Recovery functions determined from experimental data according to *Eq. 8* are often noisy because of finite sampling. Parameterized recovery functions help to overcome this problem. We tested different parameterizations (including single and double exponentials) and obtained best results with Michaelis–Menten type sigmoid functions (see e.g., Stryer 2002) (10) Here τ_{a} denotes the absolute refractory period, γ determines the maximal curvature of the recovery function, and τ_{r} is a measure of the duration of the relative refractory period in that *w*(Δ) has reached 50% of its maximum value at time Δ = τ_{a} + τ_{r}.

Parameters for the recovery function of a given cell were determined from one experiment with constant sound intensity as follows: The shortest observed ISI was taken as the absolute refractory time τ_{a}; γ and τ_{r} were then calculated from a χ^{2} fit of the theoretical ISI distribution (*Eq. 9*) to the measured data, with the stimulus strength *q* as an additional free parameter.

#### Calibration of the model.

To compare firing rates and Fano factors with model predictions, spike trains were generated by a renewal process with a bin size of *b* = 0.1 ms and conditional spike probability *b* · ρ(*t*|*t*_{last}) obtained from *Eq. 7*. The model was calibrated to the characteristics of a specific cell by computing the recovery function *w*(Δ) from one ISI distribution measured under constant-intensity stimulation at a single sound intensity *I*_{0}, as described above.

The effective input *q* of the model neuron depends on the intensity *I* of the acoustic stimulus and was determined such that the observed and the predicted mean firing rates matched. For dynamic stimuli, the procedure becomes a more elaborate 2-step process: First, the amplitude-modulated (AM) stimulus *I*(*t*) is mapped through a static nonlinearity *f*(*I*) to generate the time-varying firing rate *f*[*I*(*t*)]. Second, the effective model input *q* is calculated from *f* such that when applied to the model neuron as a constant stimulus, *q* causes the correct mean firing rate *f.*

To carry out the first step, the experimental relation *f*(*I*) is needed. Therefore firing rates *f* were measured in response to the sound intensity *I* with 10 stimuli ranging from 10 dB above *I*_{0} (the mean intensity of the AM stimulus) to 10 dB below *I*_{0}. The duration of the test stimuli was 20 ms for the highest intensities and increased linearly with decreasing intensity to 50 ms for the lowest intensities to better resolve lower firing rates. The measurement was repeated 20 times, and averages were taken (see Fig. 1). Firing-rate adaptation of locust auditory receptors causes pronounced transients in the onset response; for prolonged stimuli, the firing rate reaches a steady state after around 100–300 ms (Benda 2002). Between the test stimuli, a constant background stimulus with intensity *I*_{0} was therefore presented to keep the cell in the same adaptation state as that during measurements with the AM stimulus. The resulting “adapted” *f*–*I* curve *f*(*I*) can be parameterized by a sigmoid function derived from the positive part of a hyperbolic tangent (Benda 2002) (11) The function was fitted to the data by a least-square fit.

For the second step, the relation *q*(*f*) needs to be established. Assuming that *f* is approximately equal to the reciprocal of the median ISI, *q*(*f*) can be directly calculated for given *f*: Inserting the identity (12) into *Eq. 8* leads to *q*(*f*) · *w*(1/*f*) = 2 · *P*_{ISI}(1/*f*). Use of *Eq. 9* results in (13) Together, the model input is therefore calculated as the composition of *I* → *f* (*Eq. 11*) and *f* → *q* (*Eq. 13*), resulting in the time-dependent effective input strength *q*{*f*[*I*(*t*)]} (see also Fig. 2).

As a result of slow adaptation processes, the investigated neurons exhibited slight decreases in sensitivity over the time course of the 1-s stimulation periods. This was evident from the small but persistent decreases of the average firing rate following the strong transients characteristic for responses during the first 300 ms of stimulation. Because such long-time effects could lead to systematic errors in the analysis of short-time spike-count fluctuations, this nonstationarity was compensated as follows: The stimulus intensity used in the model was reduced by a time-dependent term such that average model firing rates in the time windows from 300 to 400 and 800 to 900 ms matched those of the recording. Because the sensitivity change was small, a linear approximation of the decay was used for simplicity. Typical corrections led to a stimulus decrease of 0.5 dB over the stimulation interval.

## RESULTS

The objective of this study is to develop a minimal model that captures essential features of spike-train variability over a wide range of input stimuli and neural response patterns, using insect auditory receptors as an easily accessible model system. Recordings were made intracellularly from the axon. We first explore the salient properties of these neurons by analyzing their responses to constant stimuli (pure tones of constant intensity). We then show that the measured spike trains are compatible with a renewal process (Cox 1962) that includes a recovery function to describe neural refractoriness (Berry and Meister 1998; Gerstner and Kistler 2002). Finally, we extend the framework to time-varying stimuli (AM pure tones) and compare the quantitative model predictions with measured test data.

Locust auditory receptor neurons encode vibrations of the tympanic membrane, the animal's eardrum, in their spike trains. In the investigated species (*L. migratoria*), firing rates depend on sound intensity in a sigmoid fashion with a maximal value at around 500 Hz (for a temperature of T = 30°C) at stimulus onset. Under prolonged stimulation, the neurons display spike-frequency adaptation. For the present analysis, we disregard the resulting initial firing-rate transient, which is most prominent during the first 300 ms (Benda 2002).

### Firing rates determine spike-train variability for constant stimuli

As a first step of our analysis, we investigated the variability of interspike intervals (ISIs) in response to sound stimuli with constant intensity. To test which parameters actually govern the ISI variability, as measured by the coefficient of variation (*CV*), we stimulated the receptors with pure tones of at least 4 different frequencies and various intensity levels. Such stimuli generally result in different firing rates because of the neuron's frequency tuning. *CV* values also vary with stimulus intensity and sound frequency (Fig. 3*A*). However, once the coefficient of variation is plotted against the observed mean firing rate, all curves coincide within the error bars, which indicates that the *CV* can be predicted from the average firing rate for constant stimuli independently of the stimulating sound frequency (Fig. 3*B*). The same observation was made for all *n* = 8 cells that were recorded under this stimulus paradigm. In general, the *CV* decreases monotonically with increasing firing rate. *CV* values near unity occur for low firing rates (<50 Hz) and approach values near 0.2 for the highest firing rates (around 300 Hz). The exact functional dependency differs from cell to cell (Fig. 3*C*).

### Spike trains are consistent with a renewal process

The intrinsic biophysical properties of a given neuron may cause arbitrarily complex spike patterns. However, if ISIs are independent under constant stimulation, the underlying neural dynamics correspond to a renewal process (Cox 1962), which allows a compact mathematical description. We therefore applied 2 particularly suited tests to investigate whether the responses to pure-tone stimuli are compatible with a renewal process (also see methods).

First, we searched for ISI correlations by calculating the serial correlation coefficients for lags of up to 20 ISIs. As shown by the data from a sample cell, the correlations are close to zero for all nonzero lags (Fig. 4*A*). The slightly negative correlation at lag 1 indicates a small effect of adaptation induced by the previous spike (for comparison also see Brandman and Nelson 2002; Chacron et al. 2001). Because the measured correlations differ by <1 SD from zero, they were neglected within the model framework.

Second, we compared measured autocorrelation functions of the observed spike train with autocorrelation functions computed from artificial spike trains with the same ISI distribution but independent ISIs (see methods). The experimentally determined and theoretically predicted autocorrelation functions closely match (Fig. 4, *B* and *C*).

Because neither test provided strong hints against assuming independent ISIs, we can proceed with this assumption and systematically explore simple renewal models to describe the measured spike-train variability.

### Recovery functions account for measured CV values and ISI distributions

All measured ISI distributions had a minimum ISI of around 1.5 ms, rose to a maximum value, and then decayed in an approximately exponential fashion (see also Fig. 5). The lack of ISIs below some minimum value and the subsequent increase reflect the influence of absolute and relative refractoriness, respectively. The exponential tail of the ISI distributions indicates that Poisson statistics is well suited for describing the generation of action potentials for sufficiently long interspike intervals.

Because the data are consistent with a renewal process, one is led to a description that accounts for neural refractoriness through a recovery function *w*(Δ = *t* − *t*_{last}). This function takes values between 0 and 1 and describes the influence of the last spike at *t*_{last} on the generation of a spike at time *t*. More precisely, the conditional rate of probabilistic spike generation ρ(*t*|*t*_{last}) is assumed to be the product of *w* (*t* − *t*_{last}) and the effective stimulus strength *q*(*t*) (14) According to this equation, spike generation is governed by 2 independent terms only, with *q*(*t*) capturing all external influences and *w*(Δ) describing a cell-intrinsic memory term. Note in particular that for constant stimuli, *q*(*t*) = *q* is a constant.

The recovery function *w*(Δ) is a unique function for each neuron and can be obtained from the cell's response to a constant stimulus, as explained in methods. We followed this approach and then used *w*(Δ) to predict the shape of ISI distributions for different sound intensities and mean firing rates through *Eq. 9*. Recovery functions were calculated from ISI distributions with an intermediate average firing rate (about 150 Hz). To avoid artifacts arising from finite sampling of the ISI distributions, recovery functions were parameterized using a class of standard sigmoid functions (see *Eq. 10*).

In the following tests of this minimal renewal model, all parameters of a cell were kept constant. To compare the model predictions to the ISI distributions derived from the test stimuli, only the stimulus strength *q* was adjusted so that the mean ISIs of the model and the recording matched. As illustrated in Fig. 5, accurate predictions of the shape of the ISI distributions can be derived with this approach. To quantify the correspondence, we compared the measured *CV* values with those calculated from the predicted ISI distributions. Figure 6*A* shows the combined data from all recorded cells (*n* = 14: 11 from the second set of experiments plus 3 with sufficient number of repetitions from the first set), each measured at 5–24 different intensities. The *CV* values from the measurements and the model match closely over the whole range of values. The variation of the observed ISI distributions and the dependence of the *CV* on the firing rate can therefore be explained in terms of the same underlying mechanism—the gradual recovery of the cell after spike generation.

The different shapes of ISI distributions from different cells can be accounted for by different parameters of the recovery functions. The steepness γ ranged from 0.4 to 5.0, with a mean of 2.4 ± 0.3. Absolute refractory periods τ_{a} were between 1.0 and 1.8 ms (mean of 1.5 ± 0.1 ms), and the parameter τ_{r} describing the relative refractory period lay between 0.8 and 4.3 ms with a mean 2.4 ± 0.2 ms. The distributions of the parameters are shown in Fig. 6*B*. Although there are 4 types of receptor neurons that differ in their attachment site to the tympanum and their frequency sensitivity (Gray 1960; Römer 1976), we did not find differences in the recovery function parameters across types. The analyzed set of cells consisted of 4 type I receptors [characteristic frequency (CF) ≅ 3.5–4 kHz], 7 type II receptors (CF ≅ 4 kHz, lower threshold than type I), one type III receptor (CF ≅ 5.5–6 kHz), and 2 type IV receptors (CF ≅ 12–20 kHz), identified by CF and best-response threshold (Römer 1976). The recovery parameters of the type I and the type II receptor neurons were not significantly different (*P* = 0.68 for γ, *P* *=* 0.64 for τ_{r}, *P* *=* 0.17 for τ_{a}; 2-sample *t*-test), and the parameters for the single type III (γ = 3.1, τ_{r} = 2.7, τ_{a} = 1.5) and the 2 type IV neurons (mean values: γ = 2.8, τ_{r} = 2.6, τ_{a} = 1.5) were also very similar to the other receptor classes.

The model developed so far provides a simple and compact description of the variability of ISIs in response to constant stimuli. As a next step, we analyzed the variability in response to time-varying stimuli and tested whether the model framework also captures response variability for this larger class of more complex stimuli.

### Recovery functions computed from constant stimuli explain the variability of responses to dynamic stimuli

Six cells were stimulated with AM pure tones at their CFs. The envelope of the stimulus was Gaussian white noise with a cutoff frequency of 400 Hz and SD of 3 or 5 dB. The mean intensity of the stimulus was set to a value *I*_{0}, which, if presented as a constant stimulus, led to an intermediate firing rate (about 150 Hz). The sound intensities thus cover the steepest region of the neuron's rate–intensity function so that the neuron is most sensitive to amplitude modulations. Matching the resulting firing-rate fluctuations as well as the ISI variability therefore presents a demanding test for the model framework.

Figure 7 shows responses of a receptor neuron that has been stimulated with such a temporally modulated acoustic stimulus. As for constant stimuli, we observe that variability is anticorrelated with the firing rate. Episodes of fast firing reflected in the high peaks in Fig. 7*B* come with low variability (Fig. 7, *E*–*G*). For counting times of 10 ms, Fano factors as low as 0.05 are found. Low firing rates, on the other hand, can lead to more than 10-fold higher Fano factors for the same counting time.

To compare the results for the time course of the firing rate and the Fano factor to the recovery-function model, we need to provide the model with a temporally modulated input *q*(*t*), which is derived at every instant from the sound intensity *I*(*t*) (see Fig. 2). As explained in detail in methods, this relation between *q* and *I* is based on the relation between *I* and the firing rate *f* (obtained from previous measurements with constant stimuli; see Fig. 1) and between *f* and *q* (calculated from the model). Drawing spikes one after the other according to the rate given in *Eq. 14*, model spike trains were generated from this input (Fig. 7*D*). As shown in Fig. 7*B*, the time-dependent firing rates of the model (black lines) and experiment (gray lines) closely coincide most of the time. This indicates that our framework of stimulus encoding, although obtained from constant-intensity episodes, is applicable to describe the responses to AM stimuli.

The model also accounts for the observed spike-count variability as measured by the Fano factor. For each of the 4 different counting window lengths, the mean Fano factors of the model (m) are very close to the experimental (e) mean Fano factors (10 ms: m: 0.25, e: 0.27; 25 ms: m: 0.18, e: 0.18; 50 ms: m: 0.16, e: 0.16; 100 ms: m: 0.15, e: 0.17; average over all cells). Moreover, the temporal fluctuations of the Fano factors are also closely met. Note that both model and measurement variability are anticorrelated with the firing rate (Fig. 7, *E*–*G*). We further assessed the quality of the Fano factor prediction by computing the mean absolute deviation of the predicted from the measured Fano factors, scaled by the error of the experimental Fano factor (see methods). The resulting relative deviations of the model predictions are near unity (mean scaled prediction error for 10-ms windows: 1.7; 25 ms: 1.3; 50 ms: 1.2; 100 ms: 1.3; average over all cells), which shows that the model deviations are of the same order of magnitude as the experimental errors of the Fano factors. Discrepancies between the measured and predicted Fano factors often coincide with mismatches between the respective firing rates (see Fig. 7*B*). This is most obvious for small counting times (Fig. 7*E*); the Fano factor is overestimated when the firing rate is underestimated, and vice versa. The population data (Fig. 8) also show that the predicted Fano factors satisfactorily match the measured Fano factors. This is true for short as well as long counting windows.

When model spike trains are compared with measured spike trains, it is notable that for sharp, high peaks in the AM of the stimulus, the model tends to exhibit greater spike-timing precision than the investigated receptor neurons (see Fig. 7*D*, e.g., at 420 ms). To quantify this, we used the correlation-based measure of spike-timing reliability by Schreiber et al. (2003). Filter widths of 0.5–3 ms were used to quantify the reliability of spike timing with different temporal resolutions (see methods). The results for all cells are shown in Fig. 9. On the shortest timescales of 0.5 and 1 ms, the spike timing of model spike trains is greater than that of the experimental spike trains. Only for 0.5-ms filter width, however, is the difference larger than the SE. For all other timescales, the reliability values are in close agreement. The deviation on the shortest timescales could be explained by additional jitter in the experimental spike trains. We recorded from the axon, so that some spike jitter might be added on the way from the spike initiation zone to the recording site, especially because the discrepancies with the model appear on the submillisecond scale. This type of jitter is not captured by our model of spike generation.

### Reduced model variants fail to capture the observed spike-train variability

To evaluate the importance of the recovery functions within the model framework, we compared it with 2 simpler, reduced model variants that are widely used in the analysis of spike-train data: rate-modulated Poisson processes without and with (absolute) refractory period. The comparison was quantified for data from a sample cell. The alternative models were tuned to this cell using the same procedures as those for the full model. The absolute recovery period was estimated as the shortest experimental ISI encountered during an experiment, and the input relations *q*(*I*) were determined by matching mean firing rates as described in methods. As expected (Gabbiani and Koch 1998; Rieke et al. 1997), Fano factors obtained from the rate-modulated Poisson process fluctuate around unity (mean Fano factor for 10-ms counting windows: 1.0; mean deviation relative to the experimental error: 16.8; see also Fig. 10*A*). Deviations from the theoretical value of unity observed in the temporal fluctuations of the Fano factor are explained by the limited number of stimulus presentations. In contrast to this, spike counts for the investigated receptor neurons display far more precision (mean Fano factor for 10-ms counting windows: 0.33).

Combining the Poisson process with a fixed refractory period to account for absolute refractoriness improves the result. The variability exhibited by this second model is closer to the measured values (mean Fano factor, 10 ms: 0.63; mean deviation relative to experimental error 7.1; see also Fig. 10*B*). However, the estimated spike-train variability is still systematically too high. This shortcoming is mirrored by errors in the ISI distribution. Only the complete recovery-function model predicts the correct amount of variability (mean Fano factor, 10 ms: 0.35; mean deviation relative to experimental error: 1.3) and the true shape of the ISI distribution (Fig. 10*C*). Thus not only absolute but also relative refractoriness is required to understand the variability and reliability of the measured spike trains.

## DISCUSSION

Stochastic models of spike generation based on recovery functions have been successfully applied to sensory neurons (e.g., see Gaumond et al. 1982; Gray 1967; Johnson and Swami 1983). Recent examples include studies about the responses of retinal ganglion cells to random flicker (Berry and Meister 1998) and responses of retinal ganglion cells, LGN neurons, and V1 neurons to moving gratings (Kara et al. 2000). In these investigations recovery functions were estimated using ISI distributions from dynamic responses to time-varying stimuli, and the model stimulus strength *q*(*t*) was inferred from the measured spike trains.

As shown by our data, responses of receptor neurons may allow an even simpler framework with a clear separation between external stimulus and stimulus-independent cell dynamics (see also Brenner et al. 2002; Johnson 1996). For the investigated auditory receptors, the cell dynamics are captured by a renewal process under constant stimulation so that each neuron can be characterized by one unique recovery function *w*(*t* − *t*_{last}). This minimal description provides accurate predictions for ISI distributions caused by input intensities over the entire range of firing-rate responses. To cover responses to dynamic stimuli, the renewal process is driven by a time-dependent effective stimulus strength *q*(*t*). The model captures the spike-count variability and salient features of the fine temporal structure in response to dynamic stimuli, although the recovery functions were always calculated from ISI distributions obtained under constant stimulation.

In the present case, application of the model was facilitated by the fact that responses from primary sensory neurons were analyzed. This allowed us to calculate the effective stimulus strength *q*(*t*) directly from the applied input and to determine the recovery functions from responses to constant stimuli. For higher-order neurons, complications arising from the dynamics of intermediate processing steps and potential feedback loops might make it impossible to fully separate internal and external contributions to the response dynamics and thus limit the present framework. However, the model could still be applied to individual modules within the sensory processing network that are then appropriately combined to describe the final input–output relation. Benefits and drawbacks of such mechanistic but multicomponent frameworks—as opposed to phenomenological one-step models such as those proposed by Berry and Meister (1998) or Kara et al. (2000)—will depend on the specific system under study.

Neural integration times have been neglected in the present model as well as in many previous approaches. Instead the model is based on an instantaneous relationship between the intensity of the acoustic stimulus *I* and the effective stimulus strength *q* so that *q*(*t*) = *q*[*I*(*t*)]. This reduction is justified by the fact that the integration time of the receptor neurons is on the timescale of about 1 ms (Gollisch and Herz 2005; Prinz and Ronacher 2002), a timescale shorter than the minimal observed ISI. This integration time is also shorter than fluctuations present in the stimuli used in this study. Stimulus integration could, however, be included in a straightforward way by describing the stimulus strength through a functional *q*(*t*) = *q*[*I*(*t*), *t*] that takes the recent history of the stimulus into account.

We did not observe differences for the parameters of the recovery function for the different receptor cell types, although we cannot completely rule out the possibility that type III and type IV cells make an exception because too few cells of these types were recorded. The 4 different types of locust auditory receptor neurons differ in their attachment site to the tympanum (Gray 1960), and their characteristic frequencies (Römer 1976), but not in their cytoanatomy (Gray 1960). Moreover, the frequency preference of a receptor stems from the resonance of the tympanic membrane at the attachment site (Michelsen 1971). This suggests that the different response classes of receptors are formed by neurons of a uniform electrophysiological type. This would lead to the observed homogeneous distribution of the recovery-function parameters over all 4 receptor-cell classes.

For testing the model's validity in response to dynamic stimuli, we focused on a restricted set of artificial stimuli. These contained substantial power at relatively high frequencies of the sound-pressure envelope because the cutoff frequency was set to 400 Hz. Stimuli with lower cutoff frequency, which are also encoded with high fidelity (Machens et al. 2001), and grasshopper communication signals were excluded, because such stimuli cause strong fluctuations in the adaptation level of locust auditory receptor neurons (Benda 2002). We rather used stimuli that fluctuate on a timescale that is much faster than the relevant adaptation time constants, which are typically around 70 ms (Benda 2002). These stimuli lead to an approximately constant level of adaptation after the initial transient at stimulus onset, allowing us to use a simplified description of stimulus encoding without adaptation. In a future extension of the model, however, one could easily include a detailed adaptation model that captures cell-intrinsic currents (Benda and Herz 2003) and dynamics of the mechanical stimulus coupling (Gollisch and Herz 2004). This would enhance the predictive power for stimuli varying on multiple timescales, such as natural grasshopper communication signals, but substantially longer recording times would be required to calibrate all model parameters. The degree of realism of the presented framework could be even further enhanced by incorporating results from biophysical investigations of the spectral (Gollisch et al. 2002) and temporal integration properties (Gollisch and Herz 2005). Altogether, this might yield a biophysically motivated and simple, yet highly accurate description of stimulus encoding for this insect auditory model system.

By comparing the recovery model to the reduced variants, we have demonstrated that both an absolute and a relative refractory period are needed to reproduce the low variance of the spike count observed experimentally. Low spike-count variances facilitate signal detection, indicating that the refractoriness of the investigated receptor neurons might be helpful for discriminating signals, such as grasshopper calling songs from different males. For this task, the temporal structure of the spike trains appears to be of particular importance. It has been shown that the spike trains of locust auditory receptors contain sufficient information for discriminating calling songs on short timescales of a few milliseconds (Machens et al. 2003). The observed renewal-process characteristics of the receptor neurons could be beneficial for this discrimination task because they enhance the coding possibilities on such short timescales. The length of an ISI generated by a renewal process depends only on the stimulus, and not on the preceding ISI, thus maximizing the number of potential output signals in the neural code. This lack of memory would be especially suited for the detection of single, fast-signal elements, such as syllables or gaps, which form parts of the communication signals of grasshoppers. On the other hand, it has been demonstrated for some systems that intrinsic ISI correlations (nonrenewal behavior) can enhance information transmission (Chacron et al. 2001). Whether a renewal or a nonrenewal process is more suitable, however, might depend on the specific task that is to be solved by the neural system.

What are the implications of our findings with respect to the neural code used by the investigated receptor cells? We have demonstrated for these cells that the mean response and its fluctuations can be predicted with a model that contains a stimulus encoder based on firing rates and a simple stochastic spike generator. Because there are about 60–80 receptor neurons per ear, which can be subdivided into 4 classes with different frequency-tuning characteristics (Römer 1976), population averages could be used to achieve reliable mean responses despite significant spike-time variability on the single-cell level. On the other hand, it has been shown that even single auditory receptor neurons contain sufficient information to discriminate conspecific communication signals, with single spikes carrying significant amounts of information at high temporal resolution (Machens et al. 2003). Moreover, some stimulus classes, especially those with a large modulation depth, are encoded much better than others (Machens et al. 2001).

This fact matches our finding that the variability of both ISIs and spike count strongly depends on the stimulus. For both constant and dynamic stimuli the variability could be modeled using the same stochastic process, which indicates that there is no principal difference between the responses to constant or strongly time varying stimuli. Differences in the encoding quality may therefore simply arise from the specific usage of the neuron's dynamic range and encoding capacity by the particular stimulus. Based on the presented framework, the interplay between specific stimulus features and the neural response dynamics could now be investigated in detail and allow a tight connection between information theory and nonlinear dynamics.

## GRANTS

This study was supported by Boehringer Ingelheim Fonds to T. Gollisch and by the Deutsche Forschungsgemeinschaft through SFB 618 “Robustness, Modularity and Evolutionary Design of Living Systems.”

## Acknowledgments

We thank H. Schütze for technical support and J. Benda and M. Chacron for helpful discussions.

Present address of T. Gollisch: Department of Molecular and Cellular Biology, Harvard University, Cambridge, MA 02138.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society