INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Variability has long attracted neurophysiologists as a tool to investigate the biophysical mechanisms of sensory processing, the integrative properties of nerve cells, and the encoding schemes used in various parts of the nervous system (Baylor et al. 1979; Hecht et al. 1942; Shadlen et al. 1996; Softky and Koch 1993). Until recently, most work has focused on characterizing the response variability of nerve cells to static stimuli, in part because simple measures such as the variance of the number of spikes recorded in long time windows provide universal and effective ways to quantify variability under such conditions (Parker and Newsome 1998).

Most biologically relevant stimuli, however, are not static. Therefore, more recently, investigators have started to characterize the trial-to-trial variability of responses to time-varying, dynamic stimuli in vivo and in vitro (Bair and Koch 1996; Berry et al. 1997; de Ruyter van Steveninck et al. 1997; Mainen and Sejnowski 1995; Mechler et al. 1998; Reich et al. 1997; Stevens and Zador 1998; Warzecha et al. 1998). When temporal variations are sufficiently strong to induce locking of spikes to stimulus transients, measures such as the standard deviation in the spike occurrence times following those transients or the probability of spike occurrence within a given time window from trial to trial may be used to provide a characterization of variability (Bair and Koch 1996; Mainen and Sejnowski 1995). However, these measures are not likely to carry over to more general stimulation conditions, when locking to stimulus transients is absent or less pronounced. An alternative approach consists of extrapolating from the study of static stimuli and to use the variance in the number of spikes observed in short time windows as a measure of variability (referred to as the spike count variance) (Berry et al. 1997; de Ruyter van Steveninck et al. 1997; Warzecha and Egelhaaf 1999). Two goals of the present work are to clarify the limits of the spike count variance as a measure of short term variability, and to introduce a new measure of spike time jitter based on recent work by Victor and Purpura (1996, 1997) that should be applicable to a wide range of stimuli, independent of the integrative properties of the investigated neurons.

Eigenmannia is a weakly electric gymnotiform fish of wave type that discharges its electric organ at regular intervals 200-600 times per second. Two types of tuberous sensory afferent nerve fibers convey information about the resulting electrical environment to the brain (Scheich et al. 1973). T-type afferent fibers provide the first stage of a pathway specialized to process phase information, called the timing pathway (Heiligenberg 1991). They fire one spike per electric organ discharge (EOD) cycle, each tightly phase locked to the zero crossings of the EOD and thus signal phase modulations of the electric field. P-type afferents, on the other hand, fire at most one spike per EOD cycle with loose phase locking to the EOD and a probability that increases in direct proportion to the mean amplitude of the field. They thus convey information about amplitude changes of the electric field to higher order neurons in the brain.

While it is well-known that the timing jitter of P-receptor afferent spikes is greater than that of T-type afferents (Scheich et al. 1973), variability in the amplitude pathway has received little quantitative attention. In contrast, variability in the timing pathway has been characterized in considerable detail, revealing the high precision of neurons in encoding phase shifts of the EOD. T-type fibers are able to fire spikes with a precision of approximately 30 µs (Carr et al. 1986). This precision increases at higher stages of electrosensory processing because of the pooling and averaging of T-type activity across the body surface (Rose and Heiligenberg 1985). Here we focus on the variability of P-type afferents and show that their firing is approximately 100 times less precise. Nevertheless, our results demonstrate that the jitter in P-receptor afferent spike trains lies within the appropriate range to efficiently convey amplitude information to the electrosensory lateral line lobe, the hindbrain nucleus that forms the first central stage of the amplitude coding pathway.

Part of our results has been presented in abstract form (Kreiman et al. 1998).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Preparation and electrophysiology

Eigenmannia specimens of 12-20 cm body length were prepared for electrophysiological recordings as described by Wessel et al. (1996) and Metzner et al. (1998). Briefly, the EOD frequency was measured, and subsequently the animal was immobilized and its EOD amplitude attenuated by an intramuscular injection of Flaxedil (gallamine triethiodide, Sigma, St. Louis, MO; <5 µg/g body wt). Under local anesthesia (2% lidocaine, Western Medical Supply, Arcadia, CA), the posterior branch of the anterior lateral line nerve was exposed just rostral to the operculum. Signals from P-receptor afferents were recorded extracellularly from this nerve with glass micropipettes filled with 3 M KCl (resistance: 40-60 M), amplified with an extra/intracellular electrometer (World Precision Instruments 767, Sarasota, FL), and stored on video tape using a PCM recording adapter (Vetter 3000A, Rebersburg, PA; sampling rate: 40 kHz). They were subsequently digitized using a commercial data analysis system (Datawave, Denver, CO; sampling rate: 10 kHz/channel). A few recordings were acquired and digitized using LabView (National Instruments, Austin, TX). Data corresponding to one point in Fig. 13 (for the cutoff frequency f_c = 88 Hz, see Stimulation below for a complete description) were obtained in a previous study (Wessel et al. 1996).

Stimulation

P-receptor afferents were stimulated as described previously (Metzner et al. 1998; Wessel et al. 1996). A sinusoidal carrier signal (Exact 519, Hillsboro, OR) with a frequency matched to the EOD frequency (f_EOD) of the fish was modulated in amplitude. The main difference with earlier work was that electric field AMs were synthesized and stored digitally for playback using commercial software (Signal Engineering Design, Belmont, MA; sampling rate: 2 kHz), allowing for repeated presentations of identical stimuli. The AM and the carrier signal were gated by the same trigger signal and were therefore phase locked to each other. The stimuli were delivered via two carbon rod electrodes, one positioned either in front of the animal or in its mouth, the other behind its tail. No differences in the neuronal responses were observed between these two configurations. The mean stimulus amplitude, measured at the side fin perpendicular to the body axis, ranged from 1 to 5 mV/cm. To avoid under-driving the afferents, it was adjusted individually for each P-receptor afferent to stimulate it at 10-15 dB above threshold.

One set of stimuli consisted of random AMs (RAMs) with a flat power spectrum (white noise) up to a fixed cutoff frequency (f_c = 5, 10, 20, 40, and 60 Hz). These AMs were obtained using a modulation signal s(t) that caused a doubling of the carrier signal amplitude for s(t) = 1 V and a reduction to zero for s(t) = 1 V (see Eq. 1 of Wessel et al. 1996). The standard deviation, , of the stimulus s(t) (which can be thought of as the stimulus contrast) was varied between 10 and 30% of the mean electric field amplitude ( = 100, 150, 200, 250, 275, and 300 mV;  = 1 V corresponded to a 100% variation of the stimulus amplitude). Consequently, amplitudes varied over a range of 20 to 10 dB of the mean stimulus amplitude. A single 15-s-long stimulus was synthesized for each parameter pair (f_c, ) and was presented 10 times, drawn in pseudo-random order from a subset of all possible (f_c, ) combinations. We usually started by presenting all f_c values at a fixed contrast ( = 250 mV) or all contrasts at two cutoff frequencies (f_c = 5, 60 Hz). Further (f_c, ) combinations were tested as time permitted.

The second set of stimuli consisted of sinusoidal AMs (SAMs) at a fixed contrast ( = 250 mV) and at various temporal frequencies f_s. The values used were f_s = 0.1, 0.5, 1, 5, 7, 10, 20, 50, 100, and 125 Hz. Each stimulus was 15 s long and was presented six times in pseudo-random order. These stimuli were presented interleaved with the RAMs protocol described above.

Characterization of spike train variability

Two methods were used to quantify inter-trial spike train variability in response to repeated presentations of the same RAM stimulus. We first computed the spike count variance as a function of the mean spike count in fixed time windows of length T (see RESULTS and Fig. 5). The same RAM stimulus was presented R = 10 times and the number of spikes, n_i, occurring in a fixed time window, T, was determined for each trial i = 1, ... R. The average number of spikes occurring in that window, n (mean spike count), and its variance, _n² (spike count variance), were estimated from
Three window sizes were used (T = 10, 50, and 100 ms), and each time window was successively shifted by 5 ms to cover the entire stimulus presentation interval. For highly variable spike trains, such as those corresponding to independent Poisson-distributed spike occurrence times, the spike count variance equals the mean independent of the window T. Conversely, if the R = 10 spike trains are exactly identical, _n² = 0 in each window T. If, however, the spike trains are not exactly identical, the minimum nonzero variance may be computed by considering the discrete nature of spiking. With f lying in the interval [0;1), we assume that a fraction (1  f) of spike counts in a fixed interval of length T equals the integer n_T (where n_T is usually small) and the remaining fraction, f, contains one additional spike, so that the spike counts equal n_T + 1. It then follows that the mean spike count, n (a positive real number), is given by

and the minimal variance is

This last equation states that _n² is a quadratic function of the fraction, f, of spike counts equaling n_T + 1 in the interval T. As a function of f, the minimal variance spans a parabola between successive integer values of the mean spike count, taking its maximal value (=1/4) at f = 1/2 and its minimal value (=0) for integer spike counts (f = 0; see RESULTS and Fig. 5). Similarly, if all spike counts in T for all R repetitions are equal to n_T or n_T + 1 except for one spike count equal to n_T  1 (or n_T + 2), then the variance still follows a parabola, but translated by a factor 2/R along the vertical axis: f(1  f) + 2/R. Successive parabolas translated vertically are generated by an analogous procedure (see RESULTS and Fig. 5).

A second measure of inter-trial variability that proved more sensitive to changes in stimulus parameters (see RESULTS) was obtained by computing an average distance between spike trains obtained in response to the same RAM. The distance measure employed was introduced by Victor and Purpura (1996) based on an earlier one used to quantify the similarity of DNA sequences (Li and Graur 1991, chap. 3; Sellers 1974). Operationally, the distance between two spike trains is defined by the following procedure: the first spike train is transformed into the second one by a series of elementary steps. Each step is assigned a "cost" and the distance is obtained by adding up the cost of all elementary steps and finding the transformation sequence yielding the minimal cost. This procedure is illustrated in Fig. 1: the two spike trains to be compared are labeled 1 and 8, while spike trains 2-7 represent the sequence of elementary steps in the transformation yielding the minimal cost. Only three elementary steps are allowed: adding a spike (as in step 6 to 7), deleting a spike (as in step 1 to 2) or moving a spike to a new position (as in step 2 to 3). The first two elementary steps are assigned an arbitrary cost of 1, whereas moving a spike by t ms is assigned a cost of q · |t| for q positive. Victor and Purpura (1996, 1997) describe an algorithm for determining the minimum cost transformation sequence and derive the mathematical properties of the ensuing distance measure, d_ij(q), between two spike trains x_i and x_j. The parameter q (measured in units of 1/time) characterizes the time interval for which the occurrence of a spike in x_i is considered to be significantly different from the occurrence of a spike in x_j: if the interval separating the spikes is larger than 2/q it is less "expensive" to transform x_i into x_j by first deleting the spike in x_i and then adding it in x_j (at a cost of 2) than by translating it to its new position (at a cost of q · |t|; Fig. 1B). It is therefore straightforward to compute d_ij(q) when q is large: let n_i and n_j be the number of spikes in x_i and x_j, respectively, and the integer c_ij denote the number of coincident spikes in x_i and x_j (coincident within some discretization interval). For large q's it is always less expensive to delete and add spikes than to move them so that the distance between x_i and x_j is obtained by first deleting (n_i  c_ij) spikes in x_i and then adding (n_j  c_ij) spikes in x_j. Thus

(1)

On the other hand, if the cost of moving a spike vanishes, q = 0, each spike in x_i may be moved at zero cost to match the position of an arbitrary spike in x_j, and a cost of 1 is only endured for each additional spike to be added or deleted in x_j. Therefore

(2)

measures the difference in the number of spikes between the two spike trains. As q  0 increases, d_ij(q) increases monotonically and reaches its maximum value (given by Eq. 1) when 2/q is smaller than the minimal time interval between two noncoincident spikes in x_i and x_j. Note that if the two spike trains are perfectly coincident d_ij(q) = 0, independent of q. The distance d_ij(q) was normalized by the total number of spikes in the two spike trains

(3)

so that d_ijⁿ(0) measures the difference in spike count normalized by the total spike count and d_ijⁿ(q  ) is the fraction of noncoincident spikes relative to the total number of spikes.

View larger version (11K): [in this window] [in a new window]   Fig. 1. Computation of spike train distances. The distance between 2 spike trains was obtained as the minimum cost to convert one spike train into the 2nd one using 3 elementary steps. A: the minimum cost path transforming spike train 1 into spike train 8 is illustrated (for a fixed value of q). Each intermediate spike train 2-7 corresponds to one elementary step: moving (from 2 to 3), adding (from 6 to 7) or deleting (from 1 to 2) a single spike. The cost of each elementary step is indicated on the right. Note that the cost of moving a spike is proportional to the distance that it is moved along the time axis. B: there are 2 alternatives to go from spike train 2 to spike train 3 in A: i) delete the last spike and add a new one or ii) move the last spike to its new desired position. The latter alternative is less expensive for the particular value of q illustrated here since q · |t1| < 2 (the dashed time interval of length 2/q corresponds to the maximum displacement for which it is less expensive to move a spike).

The effective temporal jitter, t_jitter, of the spike occurrence times was defined as t_jitter = 1/q_1/2 where q_1/2 is the value of q such that d_ijⁿ(q_1/2) = 1/2. This definition is motivated by the following arguments showing that t_jitter equals the average time interval, _inter, by which spikes are moved to transform one spike train into the other one if no spikes have to be added or deleted (see Eq. 6). Thus the effective temporal jitter t_jitter is a generalization of _inter to situations where spikes might also need to be added or deleted, as we now explain. For a fixed value of q, let n, n, and n denote the number of spikes moved, deleted, and added when computing the distance between x_i and x_j. If we pool together all noncoincident spikes in x_i and x_j, n_i + n_j  2c_ij, then each one of these spikes is either moved, deleted, or created when transforming x_i into x_j so that the following equation holds

(4)

Using Eqs. 3 and 4 to express d_ijⁿ(q) directly in terms of n, n, and n we obtain
where |t_i| is the time interval by which the ith spike (out of n) is moved. Therefore when q = q_1/2, rearranging this last equation shows that the average time interval by which a spike is moved is given by

(5)

Let us assume from now on that the number of coincident spikes is negligible, c_ij = 0 (see RESULTS). If all spikes are moved to transform one spike train into the other one (n = n = 0), Eq. 5 implies that

(6)

and 1/q_1/2 is the average time interval, _inter, by which spikes are moved. If n  0 and/or n  0, then the distance by which the remaining n spikes are moved is on average smaller to compensate for the extra cost imposed by spike additions and deletions (see Eqs. 4 and 5; the expression within the parentheses in Eq. 5 will be <1). Note, however, that the total number of displaced spikes cannot be less than half the average total number of spikes
since the right hand side of Eq. 5 has to be positive. Thus t_jitter provides an appropriate measure of spike time jitter, which automatically takes into account possible spike additions or deletions.

From the responses of a P-receptor afferent to 10 repetitions of a RAM stimulus, we computed an estimate of the average normalized distance between two spike trains as a function of q,
where n_pairs = 90 (n_pairs is obtained by considering all possible pairs of trains among 10). Normalized distances were typically computed for q = 0, 0.05, 0.1, 0.25, 0.5, and 20 ms¹ (the last value corresponds to the temporal resolution, 2/q = 0.1 ms, at which spike occurrence times were digitized). According to Eqs. 1 and 2, D_n(20) measures the average fraction of noncoincident spikes, while D_n(0) measures the average difference in spike counts (normalized by the total spike count). The average effective temporal jitter, _jitter = 1/_1/2, D_n(_1/2) = 1/2 measures the average jitter of the spike occurrence times under repeated presentation of the same stimulus. The value of _1/2 was estimated to ±0.02 accuracy [i.e., _1/2 satisfied the requirement: 0.48 < D_n(_1/2) < 0.52] by the bisection method (Press et al. 1992, chapt. 9). The percentage of spikes moved, n/(2n + n + n), and the percentage of spikes added or deleted, (n + n)/(2n + n + n), were computed over 3 s of data and six stimulus repetitions (instead of the 15 s and 10 repetitions used to compute the distances) because this task was computationally very intensive. We verified in a few cases that the results were not altered significantly by this procedure. For this latter task, a total of 15 units and 140 stimulus conditions were analyzed. We checked that the distances D_n(_1/2) computed over these reduced data sets lay between 0.45 and 0.55. This was the case for 125 stimulus conditions; the other 15 conditions were not considered further.

Stimulus estimation

The accuracy of single P-receptor afferent spike trains in encoding RAMs was assessed by linearly estimating the stimulus from the recorded spike trains. This technique essentially replaces each spike in a spike train by a continuous waveform, h(t), thus yielding an estimate, s_est(t), of the stimulus, s(t) (Fig. 2A). The waveform h(t) is chosen to optimize the match between s_est(t) and s(t) and, at low firing rates, it closely resembles the mean stimulus waveform preceding a spike (Gabbiani and Koch 1996; Wessel et al. 1996). The theoretical aspects of this signal processing technique and its application to P-receptor afferent spike trains have been discussed in detail elsewhere (Gabbiani and Koch 1998; Wessel et al. 1996; see also Gabbiani and Metzner 1999 for an introduction). For each spike train x_i(t) (i = 1, ... , 10) obtained on presentation of a RAM s(t), we subtracted the mean firing rate and estimated the filter, h_i(t), that minimizes the mean square error between the stimulus and the estimated stimulus obtained by convolving h_i(t) with x_i(t) (see Fig. 2A). This filter is called a Wiener-Kolmogorov (WK) filter in the signal processing literature (e.g., Poor 1994) and plays a role analogous to the impulse response used to estimate the instantaneous firing rate of a neuron (see Fig. 1 of Gabbiani and Metzner 1999). Each estimate of the WK filter depends on the recorded spike train x_i(t) from which it is computed and is therefore indexed accordingly as h_i(t). The WK filter was computed using MATLAB M-files (The MathWorks, Natick, MA) available at the following web address: http://www.klab.caltech.edu/~gabbiani/signproc.html. We then estimated the mean square estimation error, ², by cross-validation (Fukunaga 1990): each filter h_i(t) was convolved with a spike train x_j(t) different from the one used to compute h_i(t) to avoid over fitting. This yielded an estimate _ij²
where the brackets, ·, denote time-averaging and * denotes the time convolution operation (Gabbiani and Koch 1998). An improved estimate was obtained by averaging over all possible pairs
where n_pairs = 90. The fraction of the stimulus encoded, or coding fraction, was evaluated as
where  is the standard deviation of the stimulus. In the worst possible case, when the spike train is completely uncorrelated with the stimulus, the linear estimation algorithm predicts the stimulus mean value and the root mean square error equals the stimulus standard deviation. The root mean square error is therefore always smaller than the stimulus standard deviation (  ) so that the coding fraction, , lies between 0 and 1. The coding fraction represents the fraction of the stimulus, expressed in units of , that can be reconstructed by linear filtering of the spike train.

View larger version (20K): [in this window] [in a new window]   Fig. 2. Quantification of stimulus encoding and of its robustness to spike time jittering. A: an estimate, sest(t), of the stimulus s(t) was obtained from the spike train by convolving it with a Wiener-Kolmogorov (WK) filter (see main text for details). The accuracy of stimulus encoding by the spike train was assessed by computing the mean square error (2) between the stimulus and the estimate. The brackets, ·, denote averaging over time. B: temporal jitter was introduced by adding to each spike time a random variable taken from a zero-mean gaussian distribution with standard deviation jitter. The modified spike trains are shown for increasing values of jitter (from top to bottom) on the left. From each distorted spike train, a new WK-filter and a new estimate, sest(t), of the stimulus were computed (right). Robustness was quantified by computing the rate at which the fraction of the stimulus encoded decreased with jitter (see inset to Fig. 13). A similar procedure was used when spikes were randomly added or deleted from the spike trains.

Robustness of RAM encoding to spike time jitter, and random spike additions or deletions

To investigate the effect of spike time jitter, spike failures and the occurrence of spikes unrelated to the stimulus on the encoding of RAMs by P-receptor afferents, we created synthetic spike trains from the experimental ones by randomly adding, deleting or jittering spikes (Bialek et al. 1991). The stimulus was then estimated from these synthetic spike trains, and the coding fraction was monitored as a function of the parameters determining the amount of jitter and the number of spikes added or deleted. Each one of these three types of modifications was introduced separately. In all cases, a minimum separation of 2 ms was imposed between two spikes of the modified spike trains to take into account the refractory period of the afferent fibers.

Let p_add indicate the percentage of spikes added to the experimental spike train and p_del the percentage of spikes randomly deleted. For spike time jittering, the spikes were moved from their actual occurrence times by a random distance taken from a zero-mean gaussian distribution with various standard deviations _jitter (Fig. 2B). We used _jitter = 0, 1, 3, 5, 7, 10, 15, and 30 ms; p_add = 0, 1, 5, 10, 20, and 30%; p_del = 0, 1, 5, 10, 20, and 30%.

Let (p_add), (p_del), and (_jitter) denote the coding fractions for a given value of p_add, p_del, and _jitter, respectively. The robustness of RAM encoding by P-receptor afferent spike trains was evaluated by plotting the normalized coding fraction _n(x) = (x)/(0) where x = p_add, p_del, or _jitter as a function of x (Fig. 13, inset). In most cases, the normalized coding fraction was linearly related to the distortion parameter x (see RESULTS). We therefore performed linear fits of _n as a function of x = p_add, p_del, or _jitter


where _add, _del, and _jitter are the slopes of the regression lines. The robustness was defined as the amount of distortion required to cause a 50% drop in coding fraction
The values of p_add⁵⁰, p_del⁵⁰, and _jitter⁵⁰ were obtained by linear interpolation between adjacent values of the normalized coding fraction plotted as a function of the perturbation or by extrapolation at low stimulus cutoff frequencies (see the point f_c = 5 Hz in Fig. 13).

Modeling of P-receptor afferent spike trains

Modeling of P-receptor afferent spike trains was performed in three steps. In the first step, the variability of P-receptor afferent spike trains during RAM stimulation was compared with that of standard nonleaky integrate-and-fire models with a random voltage threshold (Fig. 3A) (Gabbiani and Koch 1998; Reich et al. 1997). The properties of the model random threshold determine the variability of the resulting spike trains. The random threshold was taken from a gamma distribution with parameters n and V_th
where
Larger values of n correspond to more reliable spike trains (see RESULTS and Fig. 7) (see also Gabbiani and Koch 1998, Fig. 9.3), and the mean voltage threshold V_th determines the mean firing rate of the model. An absolute refractory period of 2 ms was inserted after each spike occurrence. The order of the gamma distribution was varied between 1 (corresponding to an exponential distribution leading to Poisson distributed spike times), 3, 5, 10, and 100 (effectively implementing the limit n  , which corresponds to a perfect integrator). The mean voltage threshold value, V_th, was fixed so as to match the mean firing rate of the model to the one of each P-receptor afferent. Ten repetitions of the same RAM used to stimulate P-receptor afferents were fed to each model, and the distances between spike trains were computed as explained above.

View larger version (18K): [in this window] [in a new window]   Fig. 3. Comparison of P-receptor afferent spike trains to integrate-and-fire models. A: the variability of experimental spike trains was compared with the variability of perfect integrate-and-fire (I&F) neurons with a random threshold. In this model, the sum of the stimulus and a constant bias term (corresponding to the spontaneous activity) is integrated, and a spike is emitted each time that the threshold (Vthresh) is reached. After each spike, a refractory period of 2 ms is imposed and a new threshold value is chosen from a gamma probability distribution. B: to model the linear transfer properties of P-receptor afferent spike trains, the AM was first linearly filtered, with a high-pass filter fitted from the responses of P-receptor afferent to sinusoidal AMs (SAMs; see Fig. 14) and then delayed. The output z(t) was clipped and injected into a perfect I&F neuron with random threshold and refractory period equal to 2 ms.

In the second step, the linear transfer properties of P-receptor afferents were characterized using a model based on an earlier one proposed by Nelson et al. (1997) for P-receptor afferents of Apteronotus leptorhynchus (see Fig. 3B). An alternative biophysical model proposed by Kashimori et al. (1996) was not considered here, as our goal was to obtain the simplest possible description of P-receptor afferent spike trains taking into account their linear transfer function and the variability of their spike trains. The stimulus was passed through a first-order high-pass filter with transfer function H(s)

(7)

simulating the linear transfer properties of P-receptor afferents. In this equation, G_a and G_c are gain and offset terms, respectively, _a is the time constant of the filter and s = i = 2if is the complex circular frequency of the input signal. The parameters G_a, G_c, and _a were obtained by fitting the gain G(f) = |H(2if)| and the phase (f) = tan¹[ImH(2if)/ReH(2if)] of the model to experimentally measured gains and phases obtained from responses to SAMs. For each SAM stimulus, the mean instantaneous firing rate was computed over the full stimulus cycle and fitted to the function

(8)

(see RESULTS and Fig. 14). The fit parameters G_fs and _fs are the experimental gain and phase at the frequencies f_s used in the SAMs protocols, respectively (see Stimulation above). The constant c represents an offset between stimulus and response.

In the third and last step, the variability characterized in the first step and the linear filtering properties obtained in the second step were combined to obtain a complete model reproducing both the variability of P-receptor afferents and their linear filtering properties. The high-pass filtered signal was delayed by 2.5 ms (corresponding to the synaptic delay between tuberous receptors and afferent fibers), and a mean spontaneous activity was added (Nelson et al. 1997) (see Fig. 3B). The resulting signal, z(t), was then passed through a clipping nonlinearity, effectively half-wave rectifying it, and imposing a maximal firing rate of 1 spike per EOD cycle
The output, r(t) (see Fig. 3B), was fed as input to a perfect integrator with gamma-distributed threshold, as described above, to determine when a spike was fired. The order n of the gamma distribution for the threshold was selected to match the spike train variability in response to SAMs, as assessed by computing interspike interval distributions and distances between spike trains (see above and Gabbiani and Koch 1998). The responses to RAM stimuli, when available, were then compared with the model predictions (see RESULTS). In some cases the mean firing rate of the model was adjusted to take into account changes in the experimental firing rate during a recording session.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

This study is based on recordings and analysis from 69 P-receptor afferent fibers obtained in 34 different animals.

Responses of P-receptor afferents to repeated presentations of identical RAMs

To investigate the variability of P-receptor afferent spike trains and its relation to the encoding of electric field AMs, we recorded their responses to repeated presentations of identical RAMs of a sinusoidal electric field. The mean firing rates of afferent fibers were widely distributed, ranging from 25 spike/s to 374 spike/s (117 ± 69 spike/s, mean ± SD). The coefficient of variation of the interspike interval (ISI) distribution (CV = mean/SD) ranged from 0.16 to 1.7 (0.59 ± 0.36, mean ± SD). These values were similar to those observed in spontaneously active units (range: 0.12-1.12) (Wessel et al. 1996, Fig. 2B2), although several units analyzed here had higher CVs under RAM stimulation than those observed spontaneously.

Figure 4 illustrates the range of responses to repeated RAMs recorded under a variety of stimulus conditions and mean firing rates. In a few cases, the responses of P-receptor afferents were highly reproducible from trial to trial (see in particular Fig. 4C and, to a lesser extent, Fig. 4D) as has sometimes been observed in other preparations (Bair and Koch 1996; Berry and Meister 1998; Mainen and Sejnowski 1995). A clear locking of the responses to the stimulus was usually observed at high contrasts ( > 200 mV) and cutoff frequencies (f_c > 40 Hz; see Fig. 4, C and D). Furthermore, the mean firing rate of the afferent fibers had to be low (<125 spike/s; compare Fig. 4, C and G). Decreasing the cutoff frequency or the stimulus contrast tended to decrease the reproducibility of the spike occurrence times (Fig. 4, A and B). At high firing rates (>125 spike/s), P-receptor afferent responses did not show clear trends of changes in reproducibility with stimulus parameters (Fig. 4, E-H). These preliminary observations suggested that the variability across trials of P-receptor afferent spike trains depended on stimulus parameters as well as on the mean firing rate of the units.

View larger version (47K): [in this window] [in a new window]   Fig. 4. P-receptor afferent responses to RAMs exhibit a broad range of variability. A portion of the stimulus presented to each P-receptor afferent is shown on top. Each raster of spikes (9 per panel, 500 ms long) illustrates the response of the same P-receptor afferent to a single presentation of the stimulus. The left column (A, C, E, and G) illustrates responses for fixed stimulus contrast ( = 250 mV) of a neuron with low mean firing rate (A and C: mfr = 65 ± 2 spikes/s) and a different neuron with high firing rate (E and G: mfr = 137 ± 1 spikes/s) to stimuli with low and high cutoff frequencies (A and E: fc = 5 Hz; C and G: fc = 60 Hz). The right column (B, D, F, and H) illustrates the responses for a fixed cutoff frequency (fc = 60 Hz) of a neuron with low firing rate (B and D: mfr = 62 ± 1 spikes/s) and a different neuron with high firing rate (F and H: mfr = 151 ± 1 spikes/s) to stimuli with low and high contrasts (B and F:  = 100 mV; D and H:  = 275 mV).

Quantification of response variability

The spike count variance over short time windows has often been considered as an indicator of spike train variability across repeated trials of the same stimulus (Berry et al. 1997; de Ruyter van Steveninck et al. 1997). As a first step in quantifying P-receptor afferent spike train variability, we therefore plotted the spike count variance versus mean spike count across trials in windows of various sizes (10, 50, and 100 ms) as illustrated in Fig. 5. At low firing rates (Fig. 5, top row) the observed mean spike count in a given window was typically low (<10 spikes per window), and the variance across trials as a function of the mean had a scalloped appearance, reproducing almost perfectly a series of parabolas stacked onto each other along the vertical axis. Similar observations have been made in other preparations [in ganglion cells of the salamander retina (Berry and Meister 1998); in a wide-field visual tangential neuron of the fly lobula plate (de Ruyter van Steveninck et al. 1997)]. The lowest series of parabolas corresponded to the minimal possible variance that is achieved when the spike count is either equal to n or n + 1 (where n is an integer) in a given window (see METHODS). Higher parabolas corresponded successively to all spike counts equal to n or n + 1, except for one equal to n  1 (or n + 2), etc. ... per window. This result indicated that the number of spikes per window was reliable (either n, n + 1, or n  1) and well below that expected for a Poisson process (mean equals variance; dashed line in Fig. 5). However, since the scalloping was observed independently of the stimulus cutoff frequency, it did not correlate with the reliability of spike occurrence times, as observed in spike rasters (see Fig. 4, A-D). At higher firing rates, the mean spike count reached up to 25 spikes or more per window (Fig. 5, bottom row), and the variance increased considerably, ranging from the theoretical minimum up to the mean equals variance line. On average, the variance was still below that of a Poisson process.

View larger version (39K): [in this window] [in a new window]   Fig. 5. Scalloping of the variance vs. mean spike count relation is not a predictor of spike timing variability. Plots of spike count variance vs. mean spike count in windows T of 10, 50, and 100 ms. A-D were obtained in a neuron firing at low rate (mfr = 65 ± 2 spikes/s) for fixed contrast ( = 250 mV) and various cutoff frequencies fc (as indicated on the top of each panel; A is the same experiment as in Fig. 4A). E-H were obtained in a different neuron with high firing rate (mfr = 151 ± 1 spikes/s) for the same contrast and cutoff frequency values. Note that the variance vs. mean spike count curves follow the theoretical minimum curves in A-D in spite of the fact that reliable spike timing was only observed at high fcs (see Fig. 4, A-D). At higher firing rates (E-H), scalloping is still observed in some cases but is masked by a general increase in spike count variability. The 3 clusters evident in G and H, and to a lesser extent in F, correspond to the 3 window sizes (if T varies continuously between 10 and 100 ms, no clusters are observed). In all panels, mean equal to variance is indicated by a straight dashed line.

Thus according to the experimental results plotted in Fig. 5, scalloping did not appear to be directly related to the precision of spike timing across trials. To confirm this point, we artificially modified the spike trains obtained in response to repeated presentations of identical RAMs to alter the precision in spike timing without changing the statistical properties of the spike trains. We took the 10 rasters of units exhibiting scalloping of the spike count variance versus mean spike count relation and firing with varying degrees of reliability in response to RAMs (such as the rasters for the unit illustrated in Figs. 4, A and C, and 5, A-D) and successively shifted the spikes with a fixed delay t_shift. In other words, if x₁(t), ... , x₁₀(t) represent the original spike trains, new spike trains were defined as ₁(t) = x₁(t), ₂(t) = x₂(t + t_shift), ... , ₁₀(t) = x₁₀(t + 9 · t_shift). The parameter t_shift took three values: 1, 5, and 10 ms. We then computed the variance versus mean relations exactly as in Fig. 5. In all cases (5 units, 14 conditions) and irrespective of whether the timing of spikes was reliable or not, the scalloping remained present, independently of the value of t_shift. In some cases the number of vertical rows of parabolas increased with t_shift. These points are illustrated in Fig. 6, A and B. Similar results were obtained in integrate-and-fire neuron models as illustrated in Fig. 6, C and D. Thus in the worst case, t_shift = 10 ms, the timing of spikes drifted by 90 ms between the first spike train ₁(t) and the last spike train ₁₀(t) without affecting the scalloping in windows of 10, 50, and 100 ms. Since it was possible to largely eliminate any precision in the spike occurrence times from trial to trial without altering the scalloping of the spike count variance, this analysis confirmed that scalloping in these time windows was not related to the reliability of spike occurrence times.

View larger version (50K): [in this window] [in a new window]   Fig. 6. Scalloping of the variance vs. mean spike count relation measured across trials is preserved even after large shifts in the timing of individual spike trains. A: the top 10 rasters represent the response of a P-receptor afferent (mfr = 65 ± 2 spikes/s; same experiment as in Figs. 4A and 5A) to repeated presentations of a random electric field AM (RAM) stimulus ( = 250 mV, fc = 10 Hz). The corresponding spike count variance vs. mean spike count plot is scalloped as illustrated below. B: the spike trains were successively shifted by 10 ms as illustrated on top (see main text), and the variance vs. mean spike count relations was recomputed. Note that the scalloping remained present, although the variance increased as compared with A. C and D: same stimulation and analysis procedure as in A and B for an I&F neuron model with gamma order 10 (mfr = 81 spikes/s; see main text and Fig. 7A for a more detailed description of the model).

Because the spike count variance as a function of mean spike count did not offer a reliable indication of spike train variability under our experimental conditions, we turned to a second measure based on the calculation of distances between spike trains obtained under repeated RAM stimulation. This measure, D_n(q), depends on a parameter q (in units of 1/time), which determines the temporal precision at which the distance between two spike trains is computed (higher values of q correspond to higher temporal precisions, see METHODS). For two identical spike trains D_n(q) = 0 independent of q. The maximum, D_n(q) = 1, is obtained for large q's only if no spikes in the two spike trains occurred exactly at the same time. The value at which D_n(q) = 1/2, called _1/2, may be used to summarize spike train variability: if we set _jitter = 1/_1/2, then _jitter measures the average time by which spikes have to be moved to transform one spike train into the second one, or equivalently, the average jitter in spike timing. By definition, this jitter also takes into account differences in spike numbers between the two spike trains (i.e., the need to create or delete spikes to transform one spike train into the other; see METHODS and Fig. 1).

We computed the average distance between all pairs of spike trains obtained in response to the same RAM stimulus for our sample of 69 P-receptor afferents. The spike train distances D_n(q) were compared with those obtained from a family of gamma neuron models indexed by a parameter n controlling spike train variability (see METHODS). A value of n = 1 (gamma-1 neuron) corresponds to Poisson-distributed spike occurrence times in response to the stimulus while for large n (n > 100) the gamma model is identical to an integrate-and-fire neuron. Figure 7A illustrates in one example how the variability observed in P-receptor afferents compared with the model variability. The top 10 rasters labeled "P-unit" correspond to the response of a P-receptor afferent, while the next 10 rasters were obtained by simulating a Poisson (gamma-1) model. The P-receptor afferent spike trains are considerably more regular than those of a Poisson neuron and match quite well those of the gamma-10 model illustrated at the bottom of Fig. 7A. Accordingly, the average distance between two spike trains of this P-receptor afferent followed closely those of the gamma-10 neuron (see Fig. 7B,  and ) and were always smaller than those of a Poisson model (Fig. 7B, ). The value D_n(0) in Fig. 7B yields the average difference in spike number between two spike trains normalized by the total spike count. The small distance value, D_n(0) = 0.02, indicates that the number of spikes was very reproducible from one trial to the next with an average variability of 2%. On the other hand, D_n(20) is the fraction of noncoincident spikes in two spike trains at 0.1-ms resolution. The value D_n(20) = 0.98 in this experiment indicates that <2% of spikes occurred at the same time (±0.05 ms) and thus the spike trains were clearly not reproducible at a 0.1-ms resolution. The average temporal jitter in spike occurrence times, _jitter, was in this case equal to 2.9 ms (with 86% of spikes moved and 14% of spikes added or deleted), corresponding to 1.3 EOD cycles (f_EOD = 438 Hz). Furthermore, the largest deviation between D_n(q) in the gamma-10 mo