Abstract
The response of a spiking neuron to a stimulus is often characterized by its time-varying firing rate, estimated from a histogram of spike times. If the cell’s firing probability in each small time interval depends only on this firing rate, one predicts a highly variable response to repeated trials, whereas many neurons show much greater fidelity. Furthermore, the neuronal membrane is refractory immediately after a spike, so that the firing probability depends not only on the stimulus but also on the preceding spike train. To connect these observations, we investigated the relationship between the refractory period of a neuron and its firing precision. The light response of retinal ganglion cells was modeled as probabilistic firing combined with a refractory period: the instantaneous firing rate is the product of a “free firing rate,” which depends only on the stimulus, and a “recovery function,” which depends only on the time since the last spike. This recovery function vanishes for an absolute refractory period and then gradually increases to unity. In simulations, longer refractory periods were found to make the response more reproducible, eventually matching the precision of measured spike trains. Refractoriness, although often thought to limit the performance of neurons, may in fact benefit neuronal reliability. The underlying free firing rate derived by allowing for the refractory period often exceeded the observed firing rate by an order of magnitude and was found to convey information about the stimulus over a much wider dynamic range. Thus, the free firing rate may be the preferred variable for describing the response of a spiking neuron.
There has been considerable speculation about the code used by spiking neurons to transmit information (Ferster and Spruston, 1995; Sejnowski, 1995; Stevens and Zador, 1995). The spectrum of proposed theories ranges from the “rate code,” in which the firing rates of many neurons are averaged to obtain a reliable signal (Shadlen and Newsome, 1994), to “time codes,” in which the precise time relations of spikes from many neurons are meaningful (Abeles, 1991; Singer and Gray, 1995; Softky, 1995; Meister, 1996). A key factor in distinguishing among these theories is the temporal precision of individual action potentials. Thus, it is important both to measure this precision experimentally and to describe neuronal spike trains by a formalism consistent with such measurements.
The response of neurons to repeated identical stimuli is generally variable. To account for this trial-to-trial variability, such spike trains are often characterized by the instantaneous probability for generating an action potential. In this model of spike generation, termed the “inhomogeneous Poisson process” (for review, see Rieke et al., 1997), the firing probability at any instant depends only on the stimulus, not on the history of the spike train itself. Thus the variability of the response is determined by Poisson counting statistics; in particular, the variance in the spike count over any time interval is equal to the mean spike count (Rieke et al., 1997).
In recent work, we investigated the reliability of visual responses from retinal ganglion cells (Berry et al., 1997) and found far greater precision: the variance of the spike count in short time windows was much smaller than the mean. A study of the H1 interneuron in the fly visual system (de Ruyter van Steveninck et al., 1997) reached similar conclusions. In both cases, a rapidly varying stimulus caused the neurons to make sharp transitions between silence and nearly maximal firing. In addition, the spiking probability depended strongly on the history of spike times, as seen by the complete absence of spike intervals shorter than an absolute refractory period. When a neuron is firing near its maximum rate, this refractoriness causes the spike train to become more regular than a Poisson process with the same firing rate. Thus, we asked whether the refractory period plays an important role in setting the response precision of retinal ganglion cells.
Here, we compare experimental results with two different models of spike generation that modify the Poisson process in a simple way to include refractoriness. All the effects of spiking history are summarized by a recovery function, which describes how the neuron recovers its ability to fire after generating an action potential (Gray, 1967; Gaumond et al., 1982). Given a choice of recovery function, the measured spike train can be used to estimate a free firing rate, namely the firing rate of the neuron when it is not under the influence of refractoriness (Johnson and Swami, 1983; Miller, 1985). We studied stochastic models of spike generation using different forms of this recovery function and compared their predictions with experimentally measured light responses. This two-component formalism of the spiking process was highly successful at matching the observed response precision.
MATERIALS AND METHODS
Recording and stimulation. Experiments were performed on the isolated retina of the larval tiger salamander, superfused with oxygenated Ringer’s solution. Action potentials from retinal ganglion cells were recorded extracellularly with a multi-electrode array, and their spike times were measured relative to the beginning of each stimulus repeat (Meister et al., 1994; Smirnakis et al., 1997). We denote the spike train by {tij
} wheretij
is the time of theith spike on thejth stimulus trial. Spatially uniform white light was projected from a computer monitor onto the photoreceptor layer. The intensity I(t) was flickered by choosing a new value at random from a Gaussian distribution (mean
, SD δI) every 30 msec. The mean light level (
= 4 · 10−3 W/m2) corresponded to photopic vision, because the spectral sensitivity of ganglion cells equalled that of the red cone photoreceptor (Meister et al., 1994). Contrast C is defined here as the SD of the light intensity divided by the mean, C = δI/
. Most recordings used a contrast of C = 35% and extended over either 60 repeats of a 60 sec segment of random flicker or 100 repeats of a 20 sec segment.
The observed firing rate. For each ganglion cell, a peristimulus time histogram (PSTH) was obtained by histogramming the spike times tij
from all trials. Spike counts in the PSTH were divided by the observation time in each bin to yield the observed firing rate, namely the number of spikes produced per unit of time. Formally, the observed firing rate over the time interval [tA,tB
] is:
Equation 1where M = number of stimulus repeats, ρ(t) = ∑ijδ(t− tij
), and δ(t) = Dirac delta function.
The observed firing rate was calculated using time bins of 2 msec or 0.25 msec, depending on the application. For notational simplicity, we will denote it as r(t), keeping in mind that it is not a continuous function of time but is evaluated at a discrete set of time points. The same applies to all other quantities evaluated at discrete times.
Firing events. Clear firing events could be recognized inr(t) as a contiguous period of firing bounded by periods of complete silence (see Fig. 1). However, strong peaks in r(t) were not always separated by zero firing. To provide a consistent demarcation of firing events, we drew the boundaries of a firing event at minima v inr(t) that were significantly lower than neighboring maxima p1 andp2, such that
≥ 1.5 with 95% confidence. This method for identifying firing events has proven robust to changes in the parameters of the algorithm (Berry et al., 1997). With these boundaries defined, each spike in each trial was assigned to exactly one firing event.
After firing events were demarcated in this way, we characterized each event by two numbers: the time of the first spike and the number of spikes in the event. The distribution across trials of these two numbers reveals the reproducibility of event timing and spike count. For each firing event k, we computed the average timeTk of the first spike and its SD δTk across trials, which measures the temporal jitter of the first spike. Similarly, we computed the average numberNk of spikes in the event and its variance δNk2 across trials, which measures the precision of spike number. In trials that contained zero spikes for event k, no contribution was made toTk or, while a value of zero was included in the calculation of and δNk2.
RESULTS
The spike trains of retinal ganglion cells contain all the information the brain receives about the visual scene. The natural environment presents the retina with a rich variety of dynamic patterns of light. Because the response precision of a ganglion cell may not be the same for all stimulus patterns, it is important to select an experimental stimulus ensemble that provides a broad range of different stimuli. On the other hand, each particular stimulus must be repeated many times to measure the precision of the response. With these constraints in mind, we chose the stimulus ensemble provided by spatially uniform random flicker, which contains a wide variety of temporal intensity patterns (see Materials and Methods).
Previous work (Berry et al., 1997) has shown that under conditions of random flicker stimulation, retinal ganglion cells respond with discrete periods of firing separated by intervals of complete silence. These firing “events” are tightly locked to the stimulus. The time of the first spike in an event jitters very little from trial to trial; the number of spikes in the event is also very reproducible across trials. Such firing events can be identified clearly in several different ganglion cell types. They are prominent also under stimulation with spatially modulated “checkerboard” flicker. Similar results were found in the rabbit retina. In the present work, we build on the previous study by connecting the observed precision of firing events to the refractory properties of ganglion cells.
Firing events
Figure 1 shows firing events in the response of a salamander ganglion cell to random flicker stimulation. There were extensive periods in which no spikes were seen in 60 repeated trials, serving to clearly separate neighboring firing events. During such periods of firing, the firing rater(t) rose sharply from zero to a maximum value (∼200 Hz) and then fell sharply back to zero. In general, the firing events were brief bursts and could contain more than one spike (Fig. 1B).
Firing events in the response of a salamander retinal ganglion cell to random flicker stimulation at 35% contrast.A, The stimulus intensity in units of the mean during a 0.6 sec interval of the 60 sec stimulus repeat. B, Spike raster for 60 repeated trials of the stimulus. C, The observed firing rate r(t), computed by histogramming spike times in 2 msec bins.
To assess the timing precision of such a firing event, we measured the trial-to-trial jitter δT of the time of the first spike. This temporal jitter was small, typically ranging from ∼1 to 10 msec (Fig. 2A). For a given cell, δT varied somewhat among different firing events, and firing events containing more spikes generally showed greater timing precision (Fig. 2A). To characterize the timing precision of a cell with a single number, we computed the median temporal jitter τ over all events, which was 3.0 msec for the cell in Figure 2A.
The precision of timing and spike number in firing events from a single ganglion cell. During an 800 sec flicker segment repeated 30 times, 1663 firing events were identified.A, The temporal jitter δT of a firing event as a function of its mean spike count N(dots). The median temporal jitter τ is shown by the dashed line. B, The variance-to-mean ratio δN2/N of the spike count in a firing event as a function of its mean spike countN (dots), along with the value expected from Poisson statistics (thick line), and the Fano factor F (dashed line). The thin line shows the lower bound imposed by the spike count being an integer on individual trials. Note that the data points trace out a series of arches on and above the lower bound that are, successively, the next lowest possible variance-to-mean ratios.
The spike count in a firing event was also remarkably reproducible. Its trial-to-trial variance δN2 was generally less than 1 and often approached the arithmetic lower bound imposed by the fact that individual trials necessarily produce integer spike counts (Fig.2B). The ratio of the variance δN2 to the mean spike countN has a value of δN2/N = 1 for a Poisson process. By contrast, the observed variance-to-mean ratio fell below 1 for almost all the firing events and dropped below 0.1 for the largest events (Fig. 2B). The fact that δN2 << N indicates that ganglion cell spike trains cannot be characterized completely by an instantaneous firing rate (Berry et al., 1997). To assess the spike count precision of a cell with a single number, we computed the average variance over events and divided by the average spike count:F = 〈δN2〉/〈N〉. This quantity, also known as the Fano factor, was F = 0.25 for the ganglion cell in Figure 2B.
Stochastic models of the spike train
We seek an improved quantitative description of the statistics of neural firing. Toward this end, we will consider three stochastic models of the ganglion cell spike train: (1) an inhomogeneous Poisson process; (2) an inhomogeneous Poisson process with dead time; and (3) an inhomogeneous Poisson process with a renewal function. The first model is a popular approximation of spike statistics, whereas the other two models improve on the first by incorporating a refractory period. By comparing different forms of refractoriness, one can understand more readily how it affects neural precision. For each model, we show how the parameters are estimated to reproduce correctly the observed firing rater(t) (see Materials and Methods). Then we discuss how the model is used to generate simulated spike trains, the statistics of which can be compared with the measured responses.
The inhomogeneous Poisson process
We start by reviewing a very simple stochastic model of the spike train: the inhomogeneous Poisson process. Here, the probability of finding a spike in a small interval around time t is simply proportional to the instantaneous firing rate λ(t) at that time:
Equation 2In general, this firing rate varies in time, under control of the external stimulus.
Because earlier spikes do not affect the firing probability, all spikes are produced independently. Thus the probability of obtaining a particular spike train {t1,t2,… ,tn
} is proportional to the product of the individual firing probabilities:
Equation 3Under these conditions, it can be shown (for review, see Rieke et al., 1997) that the number of spikes n observed in a time interval [tA,tB] is distributed according to the Poisson distribution:
Equation 4where:
Equation 5is the average spike count. In particular, one finds from Equation4 that the trial-to-trial variance in the spike count over any given time interval is equal to the mean: δN2= N.
For this model to match the observed firing rater(t), one should choose λ(t) such that the expected number of spikes in each time bin matches the observed number of spikes. It follows from Equations 1 and 5 that this choice is simply:
Equation 6Given the firing rate λ(t), one can use this model to generate a sequence of spikes. If there is a spike at timeti
, then the probability of observing no spikes until time ti+1 is:
as can be derived from Equation 2. Thus, the next spike time ti+1 is found by selecting a random number αi+1 uniformly dis- tributed between 0 and 1, and numerically solving forti+1 the equation:
Equation 7To provide a small integration time step in Equation 7, the firing rate was sampled every 0.25 msec.
Including an absolute refractory period
It is well known that the firing probability of a neuron does depend on the history of previous spikes. In particular, there is little firing during a short period immediately after an action potential. This refractoriness results from the intrinsic dynamics of the active membrane conductances that are responsible for spike generation (Hodgkin and Huxley, 1952).
The simplest method of incorporating a refractory period into a stochastic model of spike generation (Gray, 1967; Gaumond et al., 1982;Johnson and Swami, 1983) assumes that the firing rate λ(t,{tij
}) of a neuron is the product of a free firing rate q(t) and a recovery function w(t −tlast) that reduces the firing rate immediately after a spike:
Equation 8The free firing rate q(t) is a function of the stimulus, whereas the recovery functionw(t − tlast) depends only on the time since the last action potential attlast. It varies from 0 at short times to 1 at long times.
We begin by considering an absolute refractory period: a neuron that has fired an action potential is unable to produce another for a period of time μ, regardless of the strength of the stimulus. After this dead time, its firing rate returns immediately to the value in absence of refractoriness. Formally, the recovery function is:
Equation 9
where
is the Heaviside step function.
Given a choice of dead time, we would like to estimateq(t) from the spike trains so that this spiking mechanism reproduces the observed firing rater(t). At time t during the stimulus repeat, the model cell fires at the instantaneous rateq(t), but only if t does not fall into the dead time of the previous spike. The probability that this happens can be estimated directly from the measured spike trains, as illustrated in Figure 3. On a single trial j, define:
Equation 10which can be expressed formally as:
Equation 11Then the probability that time t is available for free firing is given by the average ofWj
(t) over all trials:
Equation 12The observed firing rate of the model cell equals the free firing rate q(t) multiplied by the probabilityW(t) that free firing is possible at timet. Equating this with the observed firing rater(t) leads to:
Equation 13The general form of Equation 13 was introduced by Johnson and Swami (1983), and was later shown by Miller (1985) to be the maximum likelihood estimate for q(t). An alternative method estimates the free firing rate directly from the PSTH using an iterative procedure (Jones et al., 1985; Bi, 1989).
Schematic illustration of the derivation ofW(t), the probability of free firing. On each trial j,Wj (t) (solid gray) has a value of zero during the refractory period following a spike (vertical lines) and one elsewhere. This function is averaged over all trials (bottom), to yield W(t), the fraction of trials during which free firing was possible at time t.
The fraction W(t) depends on the choice of refractory period, as seen in Figure 3. Thus, the estimate of the free firing rate will be different for each choice of μ. To make this dependence explicit, we denote the free firing rate byqμ(t). In particular, when there is no refractory period, then W(t) = 1. Thus the nonrefractory model defined by Equation 6 is seen as a special case of Equation 13: λ(t)=qμ=0(t). For a non-zero dead time μ, W(t) ≤ 1, and the free firing rate obeys the inequality q(t) ≥r(t). If the refractory period μ is chosen too large, it may exceed some of the interspike intervals in the data. This can lead to a complication in estimatingqμ(t): if a spike was observed at time t, but this time lies within μ of the preceding spike on every trial, then W(t) = 0 andqμ(t) = ∞. To avoid the divergence, we imposed a maximum bound on the free firing rateqμmax(t) = 1000r(t). As described below, models in which the chosen refractory period μ is too large eventually fail to account for the observed firing rate r(t).
Given the free firing rate q(t) and the recovery function w(t), one can again generate simulated spike trains from a set of random numbers {αi}, uniformly distributed between 0 and 1. If there is a spike at timeti, then the next spike timeti+1 is found by numerically integrating:
Equation 14which follows by substituting the firing rate λ(t,{tij
})from Equation 8 into the spike-generating rule of Equation 7. Again, for bothq(t) and w(t) we used an integration time step of 0.25 msec. For an absolute refractory period, Equation 14 simplifies to:
Equation 15
Including a relative refractory period
In practice, the firing probability of a neuron does not recover instantaneously after the dead time has expired. In a more realistic model of spike generation, each action potential is followed by a recovery period, during which the firing rate is gradually restored to its free value (Kuffler et al., 1957). This amounts to choosing a recovery function w(t) that approaches 1 gradually, rather than in the step-like manner considered above. A good choice for w(t) can be derived directly from the observed spike trains, as follows.
Assume that the neuron behaves according to Equation 8 and that its free firing rate is constant, q(t) =q. Then the time interval Δ between successive spikes is distributed according to:
Equation 16For long interspike intervals, where w(Δ) ≈ 1, one predicts an exponential dependence p(Δ) ∝ exp(−qΔ), whereas the frequency of short intervals, wherew(Δ) < 1, should fall below this exponential curve. The first term in Equation 16 is the probability that there is no spike in a time interval Δ, equal to: 1 − ∫0Δ p (Δ′)d Δ′. Thus, one can express the recovery function in terms of the interspike-interval distribution:
Equation 17Figure 4A shows the distribution p(Δ) of interspike intervals for a ganglion cell under random flicker. For long intervals between 5 and 10 msec, the distribution drops exponentially, as expected for a Poisson process with a constant firing rate. The complete absence of very short intervals indicates that absolute refractoriness lasts for 2 msec. The relative paucity of intervals between 2 and 5 msec results from a relative refractory period, during which the ganglion cell is able to fire an action potential, but does so with reduced likelihood. To determine the recovery function w(t) that governs this period, we first estimated q as the rate of the exponential decay at long times (Fig. 4A) and then calculated w(t) for intermediate times from Equation 17. This function is shown in Figure 4B. Note that w(t) is negligible out to 2.5 msec and then rises monotonically between 2.5 and 5 msec.
Illustration of the method for determining the recovery function wm (t) from the spike train of a ganglion cell. A, The histogram of interspike intervals Δ (diamonds) is fit over the range Δ = 5–10 msec with an exponential curve P(Δ) ∝ e−qΔ(solid line), yielding a decay rate ofq = 780 Hz. B, This value ofq serves to compute the recovery functionwm (t) using Equation 17. For intervals >10 msec, P(Δ) was approximated by the exponential extrapolation shown in A.
The derivation of Equation 17 assumed a constant free firing rate, whereas the observed firing rate was clearly not constant, instead showing rapid modulations (Fig. 1). However, the shape of the interspike-interval distribution at short times is dominated by the rapid firing during firing events. The peak firing rates observed did not vary much from event to event, typically within a factor of 2 for events having at least one spike per trial. Consequently, deriving the recovery function in this way may still provide a good estimate of how spiking recovers. This was confirmed by the practical success of this model, as described below.
Given this choice of w(t), we derived the free firing rate by the general procedure described above. First, the effective probability of free firing W(t) was computed from w(t) and the spike trains (Eq. 11and 12). Then the observed firing rate was divided by the probability of free firing to yield the free firing rate q(t) (Eq. 13). Because w(t) is derived directly from the interspike-interval distribution, its absolute refractory period extends only over a duration for which there are no interspike intervals in the observed spike trains; thus these calculations were never hampered by the divergence of q(t) discussed earlier. To identify the recovery function and free firing rate associated with this stochastic model incorporating the most realistic refractory period, we will denote them aswm (t) andqm (t), respectively. This model again produces simulated spike trains from random numbers by the general procedure of Equation 14.
Performance of the models
In the following section, we will test the predictions of the models developed above against the observed light response of retinal ganglion cells. Simulated spike trains were produced with the same number of repeated stimulus trials as during experiments and analyzed in the same manner (see Materials and Methods). The firing rater(t) was computed from the PSTH, firing events were identified, and both the temporal jitter of the first spike of an event δT and the variance in the spike count during the event δN2 were calculated. To further test the models, four other statistics of the simulated spike trains, described below, were compared with the corresponding values for the observed spike trains. Quantities derived from simulations will be denoted by subscripts, for examplerμ(t) for the firing rate of the model with absolute refractory period of length μ, andrm (t) for the model with a relative refractory period.
Spike train statistics
Two statistics were used to assess the accuracy of the firing rate produced by the model rμ(t). The simulated firing rate averaged over the entire stimulus segment
μ was compared with the observed mean firing rate
. To test how well the model captured the fast modulations in the observed rate, we calculated the mean-squared error of the predicted rate normalized to the total variance in the observed rate:
Equation 18For this purpose, r(t) was evaluated with a bin size of 2 msec, less than the temporal jitter τ for a typical ganglion cell, and thus fine enough to capture the steep rising and falling edges of firing events. Because the responser(t) varies somewhat from trial to trial, it is of interest to see how the quality of the model’s fit compares with this intrinsic variability. The corresponding benchmark of the response variability is given by:
Equation 19where δr(t)/
= standard error of r(t) across trials.
Similarly, the random error associated with the model’s prediction is measured as Eμ0, obtained by replacing rμ for r in Equation 19. In particular, when Eμ =Eμ0, then the model agrees with the observed firing rate to within the uncertainty in the model’s firing rate, which is limited by the finite number of trials.
The event statistics δT and δN2 both measure the trial-to-trial reproducibility of spike trains, but only in very specific aspects. Thus, we also computed a much more general measure of variability, the so-called “noise entropy” of the spike trains (de Ruyter van Steveninck et al., 1997; Strong et al., 1997). In brief, the spike train was evaluated in 2 msec bins, and thus converted into a string of ones and zeroes marking the presence or absence of a spike. This string was evaluated one “word” at a time, with word lengths ranging up to 25 bins. At a given time t during the stimulus segment, theM trials yielded M words, each providing the local spike pattern during one trial. The variability of these spike patterns was measured by the entropy of the set of words (Shannon and Weaver, 1963). By averaging this measure over all times t, and dividing by the duration of a word, one obtains the rate of noise entropy per unit time, Snoise. Using 60 or 100 trials (see Materials and Methods), the variety of words could be sampled well for words up to 50 msec long. As a result, this measure captures the variety of spike patterns within a firing event.
The noise entropy measures how much the local spike pattern varies from trial to trial. It is important to compare this variability introduced by neural noise to the overall variation in spiking patterns elicited by the stimulus. Therefore, we also computed the total entropy of the spike train, Stotal. This is derived in a similar manner from the set of all words encountered across times throughout the stimulus segment. Finally, the difference between the total entropy and the noise entropy, I =Stotal −Snoise, represents the information conveyed by the spike train about the visual stimulus (Strong et al., 1997). This can be seen as follows: without knowledge of the stimulus, the uncertainty about what spike pattern the cell might produce in the next 50 msec is equal to the total entropyStotal. If the stimulus is known, this uncertainty is reduced to the noise entropySnoise. The difference is equal to the information gained about the spike train given the stimulus, which by symmetry is equal to the information about the stimulus gained from the spike train.
Predictions from an absolute refractory period
To gain a basic understanding of how refractoriness affects response precision, we first analyzed the effects of an absolute refractory period. We varied the absolute refractory period μ from 0 to 6 ms, and analyzed how the simulated spike trains produced by the model changed. Figure 5 shows the results for a single representative retinal ganglion cell. The simulated average firing rate
μ matched the actual firing rate of the neuron very well up to refractory periods of μ ≈ 4 msec (Fig. 5A). For larger values of the absolute refractory period, many interspike intervals in the responses of the cell are shorter than μ. In those cases, the model cannot match the firing rate of the cell, and thus the predicted rate drops, although the discrepancy at μ = 6 msec is still only 10%. Figure5B compares the mean-squared errorEμ of the model with the uncertaintyEμ0 in its predicted firing rate. For refractory periods in the range to msec, the mean-squared errorEμ comes very close toEμ0 and also to the uncertaintyE0 in the observed firing rate of the neuron. This suggests that the model matches the time course ofr(t) about as well as can be expected from the finite number of stimulus trials; no systematic deviations from the measured firing rate can be resolved in this regime. Above μ = 4 msec, the mean-squared error rises, because the predicted firing rate is systematically too low, as discussed above. Below μ = 1 msec, the mean-squared error Eμ also rises, although it still agrees with the value Eμ0 expected from the trial-to-trial variation. This occurs because the variability of a nonrefractory spike train (μ = 0) is greater than that of the real neuron.
Comparison of the model with an absolute refractory period to observed results from a single representative ganglion cell. Six statistics of the spike trains produced by the model are plotted as a function of the absolute refractory period μ. The value of each statistic was averaged over 10 sets of simulated spike trains, and error bars denote ± 1 SD. In each panel, the real value from the neuron is shown by a dashed line.A, The average firing rate
μ
; B, the mean-squared error Eμ (solid symbols) and the uncertainty Eμ0(open symbols) in the simulated firing rate, with the corresponding uncertainty E0 in the measured rate (dashed line); C, the Fano factor Fμ; D, the temporal jitter τμ; E, the noise entropy per unit time Sμnoise; F, the information per spike Sμtotal−Sμnoise)/
μ.
Although the mean firing rate was approximately constant for refractory periods up to 4 msec, the precision changed dramatically. Figure5C shows that the Fano factor Fμvaries from the expected value of 1 for no refractory period to below 0.2 for the largest refractory period. This happens because refractoriness leads to a regular spacing of spikes during periods of rapid firing, and this reduces the variability in the number of spikes generated during an event. The spike number precision of the model matched the observed value of F = 0.25 at a refractory period of μ = 4.5 msec. The temporal jitter τμ also decreased as refractoriness was added (Fig. 5D), although the effect was somewhat smaller. This sharpening of temporal precision occurs because the free firing rateqμ(t) rises more steeply thanr(t) (see Fig. 8), so that the first spike of an event occurs over a narrower range of times. The model’s timing precision matched that of the neuron at μ = 2.5 msec and showed a 10% discrepancy at μ = 4 msec.
The improved response precision from an absolute refractory period resulted in higher information transmission rates. The noise entropySnoise remained constant until μ = 2 msec, but then decreased sharply (Fig. 5E). Thus, refractoriness significantly reduced the variety of different spike patterns produced during a given firing event. The total entropyStotal also decreased, but by a smaller amount, because most of the variety of spiking patterns across events is related, in fact, to variety in the stimulus. Consequently, the visual information rate per spike (Sμtotal−Sμnoise)/
μincreased by ∼30% as μ ranged from 0 to 6 msec (Fig.5F). Both the noise entropy and the information rate of the model matched the real values for refractory periods of μ = 4 msec, and the same held to good approximation for the other comparisons tested here. Therefore, a probabilistic spike generator with an absolute refractory period can reproduce many statistics of these ganglion cell spike trains.
Predictions from a relative refractory period
The most evolved model that incorporates a relative refractory period with a smooth recovery functionwm (t) was tested in a similar manner, but the results will be elaborated in greater detail. For each cell,wm (t) andqm (t) were determined from the responses as described above; recall that this entails no free parameters. Then the statistics of simulated spike trains were compared with those of the actual response.
Figure 6 compares the precision of firing events simulated by the model with those in the real responses of a single ganglion cell. The temporal jitter δTm matches the real δT rather well. There may be a barely detectable downward bias for the simulated jitter δTm . Similarly, the variance δN2 in the spike count during events is matched by the model’s prediction δNm2over a wide range of values.
Comparison of the model with a relative refractory period to observed results from a single representative ganglion cell. A, The predicted temporal jitter δTm plotted as a function of the observed value δT for 213 firing events (dots).B, The predicted spike number variance δNm2 plotted as a function of the observed value δN2 (dots). In each panel, the identity relation is shown as a dashed line.
To better assess the validity of this model, we studied its performance on a population of 42 salamander ganglion cells in three retinae. The average firing rate produced by the model was very close to the observed value: the discrepancy between predicted and real firing rates ‖
m
−
‖/
was 1.6% on average over the population of cells. The mean-squared error Em
of the predicted firing rate modulation again came very close to the intrinsic uncertaintyE0 in the response of the cell: the average ratio wasEm
/E0 = 1.1. Thus, the model captures the time course of the firing rate as well as possible, given the finite number of trials.
For each cell, the precision of firing events was again measured by the Fano factor F and the temporal jitter τ. Figure7A demonstrates that the Fano factor Fm predicted from the model agrees very well with the observed number precision F for almost all cells. By contrast, the nonrefractory spike generator always produces a Fano factor Fμ=0 near unity, as expected, and in stark disagreement with the real neurons. Figure 7B shows a similar comparison for the temporal jitter of firing events. Here, one finds that both a relative refractory period and the simple nonrefractory model show close agreement with the measured temporal jitter, with a slight negative bias using the relative refractory period and a slight positive bias for the nonrefractory model. This similarity is expected, because the timing precision is derived from the first spike of an event, which is least affected by refractoriness.
Performance of the model with a relative refractory period and the nonrefractory model for a population of 42 ganglion cells from three retinae. In each panel, a statistic derived from simulated spike trains is plotted against the observed value from the same cell; the identity relation is shown as a dashed line. A, The Fano factorFm (open circles) andFμ=0 (solid squares) as a function of F. B, The temporal jitter τm (open circles) and τμ=0 (solid squares) as a function of τ. C, The discrepancy in the noise entropy ΔSmnoise =Smnoise −Snoise (open circles) and ΔSμ=0noise= Sμ=0noise −Snoise (solid squares) as a function of Snoise.
The model with a relative refractory period also matched the variety of spiking patterns as measured by the spike train entropy. The total entropy Smtotal predicted using a relative refractory period was similar to the real valueStotal: the discrepancy was ‖Smtotal −Stotal‖/Stotal= 2.9% averaged over the population. Similarly, the noise entropy was predicted very well with a relative refractory period (Fig.7C). By contrast, the nonrefractory model produced a large discrepancy between predicted and observed noise entropy, ranging up to 5 bits/sec. In conclusion, a spike generator that incorporates a relative refractory period can very closely reproduce many statistics of real spike trains, whereas a nonrefractory model clearly fails.
The free firing rate
In this model of the spike train, the free firing rate depends only on the external stimulus, whereas the observed firing rate is separated from the stimulus by an additional, history-dependent nonlinearity. As a result, the free firing rate may be a more fundamental response measure than the observed firing rate. In addition, a refractory period causes the observed rate to saturate whereas the free rate is not so constrained. It is therefore conceivable that the free firing rate can distinguish gradations in the visual stimulus that are lost in the observed rate. Motivated by these ideas, we studied the behavior of the free firing rate of retinal ganglion cells under random flicker stimulation. The free firing rate was estimated using the model with a relative refractory period, which emerged from the analysis presented above as the most realistic implementation of the refractory properties of a ganglion cell.
Comparison with the observed firing rate
Figure 8 compares the free firing rate qm (t) with the observed firing rate r(t) for a typical firing event. The observed rate rises sharply from zero to its maximum of ∼300 Hz and then exhibits several peaks (Fig. 8A). The simulated firing rate rm (t) closely follows these peaks. Figure 8B illustrates the free firing rate qm (t) on the same scale. At the beginning of a firing event, when refractoriness is negligible,qm (t) is equal tor(t), but it becomes much larger after several spikes have occurred (Gaumond et al., 1983), ranging up to ∼3000 Hz in this case. As a result, the waveform of is even narrower than that of r(t). Furthermore,qm (t) is smoother thanr(t), because the ripples introduced by refractoriness are largely removed. Again, this suggests thatqm (t) is more closely related to the excitatory input that this neuron receives.
A, The observed firing rater(t) (solid line) and the simulated firing raterm (t) (dashed line) generated by the model with relative refractory period for a single firing event. B, The observed firing rater(t) (thin line) and the free firing rateqm (t) (thick line) for the same firing event. Each curve was computed with 0.25 msec bins, for greater time resolution, and then boxcar-smoothed over nine bins to roughly correspond to the 2 msec bins used elsewhere. C, The observed firing rater(t) computed in 2 msec bins plotted against the corresponding values of the free firing rateqm (t) (dots). Also shown is the steady-state relation r =q/(1 + qμ) for absolute refractory periods of μ = 1.5 msec (dashed line) and μ = 3.5 msec (solid line).
Figure 8C is a scatter plot of the values ofr(t) against values ofqm (t) at the corresponding time throughout the entire stimulus period. Notice that they are roughly equal at low firing rates (q ∼ 10 Hz), but that r saturates at its maximum of ∼300 Hz whereas q still increases. The general shape of this relationship can be understood as a simple consequence of refractoriness: if q(t) were constant and there were an absolute refractory period of μ, then r would reach a steady-state value of r = q/(1 +qμ). This function provides an upper envelope to the scatter plot in Figure 8C with μ chosen as 1.5 msec, close to the absolute refractory period of the real data. Matching the middle of the data cloud is a similar curve using μ = 3.5 msec, roughly the value that best reproduced the spike train statistics.
The free firing rate qm (t) still varies by more than an order of magnitude under conditions where the observed rate r(t) has already saturated (Fig.8C). However, one may ask whether this variation is systematic or is instead introduced by noise in the algorithm that estimates qm (t). After all, the free firing rate is given by an expression (Eq. 13) with a denominator that is close to zero for large values of q. To examine this concern, we tested whether these variations in the free firing rate are driven systematically by the light stimulus.
The light response of the free firing rate
For concreteness, we will assume that the ganglion cell firing rate r(t) depends on the stimulus intensity with what is called a “Wiener cascade” or “LN model” (Hunter and Korenberg, 1986):
Equation 20where
Equation 21Here, the stimulus I(t) is first convolved with a linear filter L(t), and the resultZ(t) is transformed at each point in time by a nonlinear function f(Z). This simple LN model has been used with some success in previous analyses of visual responses (Korenberg et al., 1989; Mancini et al., 1990; Sakai et al., 1995). For a mechanistic picture, one could interpretL(t) as a temporal filter implemented by retinal circuitry. The filtered stimulus Z(t) can be thought of as the effective input to the ganglion cell, and the function f(Z) transforms this input into a firing rate. Of course, the following analysis does not depend on this interpretation, nor does it claim that the LN description is a complete account of the light response of the ganglion cells.
The parameters of the LN model were deduced from the ganglion cell spike train as follows. First, the linear filterL(t) was obtained by reverse correlation ofr(t) to the Gaussian stimulusI(t) (Hunter and Korenberg, 1986):
Equation 22Then, the effective input functionZ(t) was computed from Equation 21. Finally, the transform f(Z) was found by plotting the measured firing rate r(t) against the effective input Z(t). The same analysis was applied to derive the LN relationship that links the free firing rateqm
(t) to the stimulus. Obviously, this yielded different parameters, which will be denoted asLm
(t) andfm
(Z).
Figure 9A shows an example of the two linear filters derived from the observed firing rate and the free rate. The large downward component preceding the action potential indicates that this cell was OFF-type. Both filters have a standard biphasic shape (Sakai et al., 1988; Meister et al., 1994;Smirnakis et al., 1997), implying a band-pass frequency characteristic, and thus a transient light response. The filter from the free firing rate Lm (t) is somewhat stronger thanL(t), especially in its second peak, and also has slightly shorter latency. This enhancement of the filter can be explained if the largest values of q follow more intense stimulus patterns that involve larger excursions from the mean light level. Stronger stimuli also tend to elicit responses faster, consistent with the shorter time-to-peak forLm (t).
Analysis of the light response for the observed firing rate r(t) and the free firing rate qm (t).A, The linear filter relating the light stimulus to the effective input: L(t) for the observed firing rate (thin line) andLm (t) for the free firing rate (thick line). B, The response function relating the effective input to the firing rate:f(Z) for the observed firing rate (thin line) andfm (Zm ) for the free firing rate (thick line).f(Z) was obtained by plotting values of r(t) againstZ(t) every 5 msec, and then binning the range of Z at intervals of ΔZ = 0.4 SD. Thus,f(Z) is the average value ofr(t) over all time points that have an effective input Z(t) between Z − ΔZ/2 andZ + ΔZ/2. The same procedure was followed with qm (t) andZm (t) to obtainfm (Zm ).
Both linear filters were convolved with the stimulus to obtain the effective input Z(t) andZm (t). Large values of the effective input occur when the preceding stimulus had a time course similar to the linear filter, whereas small values arise from stimulus patterns that are unrelated. A positive value of Z marks an intensity variation that matches the sign of the linear filter, for example a brief dimming of the light ∼70 msec earlier (Fig. 9A), whereas a negative value indicates an intensity variation of the opposite sign.
Figure 9B shows the response functionf(Z) obtained by plotting the firing rater(t) against its effective inputZ(t). At negative values of Z, the relationship produces a low floor of firing at ∼0.1 Hz, which appears unrelated to the effective input. Clearly, this cell does not convey much information about intensity transients of the “wrong” sign. At positive values of Z, the firing rate rises steeply, and follows an exponential dependence over almost three orders of magnitude. Then f(Z) saturates at a firing rate of ∼100 Hz. The corresponding response function for the free firing ratefm (Z) follows a similar dependence at low Z. However, unlikef(Z), it does not saturate at high values of Z, but continues to increase for more than another order of magnitude to ∼2000 Hz. Thus, the free firing rateqm (t) can report the strength of the stimulus over a much wider range than r(t), up to the strongest intensity fluctuations encountered in these experiments.
DISCUSSION
In summary, we have shown that the light response of a retinal ganglion cell can be captured well by its free firing rateq(t) and its recovery from refractorinessw(t). Both of these functions can be estimated directly from the raw spike trains without free parameters. This two-component description of the spiking process byq(t) and w(t) has significant advantages over the traditional analysis of the observed firing rate r(t): it accounts correctly for the trial-to-trial variation of the response, the precision of spike timing and spike numbers, and other statistics of the spike train. Furthermore, q(t) continues to vary under conditions in which r(t) is saturated, and continues to distinguish gradations in the response of the neuron. Thus, q(t) may prove useful for interpreting the response of any spiking neuron.
Our method relies on determining the recovery functionw(t) directly from the distribution of interspike intervals. The recovery function includes both an absolute and a relative refractory period that extend over a short time interval of ∼5 msec after a spike. These characteristics are consistent with the refractory properties of the ion channels that produce action potentials: sodium channel inactivation leads to an absolute refractory period and relative refractoriness can result from either shunting or hyperpolarization while the potassium channel remains open (Hille, 1984). Because these ion channels are ubiquitous in spiking neurons, one expects the present results to be widely applicable also outside the retina.
In fact, the effects of refractoriness on signaling have been studied in cochlear nerve fibers (Gray, 1967; Gaumond et al., 1982, 1983;Johnson and Swami, 1983; Miller and Mark, 1992). Their spike trains reveal absolute and relative refractory periods similar to those of salamander retinal ganglion cells. In cochlear afferents, refractoriness was found to suppress the observed firing rate relative to the free rate by a factor of 2–5 for acoustic clicks and tone burst stimuli (Gaumond et al., 1983; Miller and Mark, 1992). At the onset of spiking episodes, the observed and free firing rates agreed, but the observed rate became most strongly suppressed when it achieved its highest values (Gaumond et al., 1983). In related studies, refractory spike generation helped to explain the power spectrum of spike trains from cochlear afferents (Edwards et al., 1993) and from neurons in area MT of the visual cortex (Bair et al., 1994).
Spike trains may also show a long-term dependence on previous activity. Cochlear afferent fibers have, in addition to a refractory period lasting ∼4 msec, a gradual recovery period extending 20–40 msec after an action potential (Gray, 1967; Gaumond et al., 1982).Lowen and Teich (1992) found that successive interspike intervals were not independent, but instead showed negative correlations over time intervals up to 100 msec. When firing spontaneously or responding to sustained tones, auditory neurons exhibit a spiking variance on long time scales greater than expected from a Poisson process (Teich et al., 1990). Furthermore, salamander ganglion cells show a slowly inactivating outward membrane current (Lukasiewicz and Werblin, 1988) that could contribute to a gradual recovery from spiking. However, when driven by a visual stimulus, these neurons appear strongly locked to the intensity modulations. Spikes are grouped in well defined firing events, and we have found little evidence for correlations extending beyond the typical duration of a firing event. Thus, the present work focused on a recovery function w(t) that implemented only short-term refractoriness.
How should the free firing rate q(t) be interpreted? It is tempting to imagine that q(t) is directly related to the synaptic currents entering the ganglion cell, whereas the saturation in r(t) results from the refractory properties of its active membrane. If this were the case, then the relationship between q and robserved here (Fig. 8) could also be measured by injecting current into a ganglion cell through an intracellular electrode. Such experiments with retinae from salamander (Diamond and Copenhagen, 1995) and turtle (Baylor and Fettiplace, 1979) have reported a rather linear dependence between the injected current and the ganglion cell firing rate. A similar linear relationship has been observed in other spiking neurons (Granit et al., 1966; Lanthorn et al., 1984; Powers et al., 1992). However, the peak firing rates achieved by current injection (50 Hz in salamander) (Diamond and Copenhagen, 1995) were not high enough to sample refractoriness. Clearly, it would be rewarding to explore the dependence of firing rate on synaptic input over a greater range.
Our analysis of the effects of refractoriness indicates that a longer absolute refractory period serves to improve the response precision of ganglion cells and increase their information rates. Similarly, in the auditory system, realistic refractoriness has been shown to reduce the variance in the response of cochlear afferents to vowel sounds from its value with no refractoriness (Miller and Mark, 1992). These results suggest that a refractory period can improve neural signaling. The refractory period is commonly thought to limit the performance of a neuron, but this view is tied to a preconception about how the cell conveys its signal. If neurons use a “firing rate code,” in which the message lies in the instantaneous frequency of firing, then refractoriness does reduce the dynamic range of neural output by causing saturation in the firing rate. However, if the message lies in the number of spikes fired during a discrete firing event, then the maximum firing rate is not a limitation, because the cell can simply fire for a longer period of time to produce more spikes. In this case, the refractory period improves the temporal precision of the second and subsequent spikes in an event, leading to a spike count of much higher fidelity than would obtain from a Poisson process. Considering the receiver of such nerve messages, it is plausible that postsynaptic neurons “count” the number of spikes in an event, because the typical firing event duration of ∼20 msec is within a postsynaptic integration time, whereas the interval between events is much greater.
Footnotes
This work was supported by a National Research Service Award to M.J.B., and a Presidential Faculty Fellowship and Pew Scholar Award to M.M. We thank Mike DeWeese for many useful conversations.
Correspondence should be addressed to Michael J. Berry II, Department of Molecular and Cellular Biology, Biological Laboratories, 16 Divinity Avenue, Cambridge, MA 02138.
