Modeling the Variability of Firing Rate of Retinal Ganglion Cells MICHAEL
W. LEVINE
Uniuersityof Illinois at Chicago, Department of Psychology and Committee on Neuroscience,
Chicago, Illinois 60680 Received 4 October 1991; reuised 28 July 1992
ABSTRACT Impulse trains simulating the maintained discharges of retinal ganglion cells were generated by digital realizations of the integrate-and-fire model. If the mean rate were set by a “bias” level added to “noise,” the variability of firing would be related to the mean firing rate as an inverse square root law; the maintained discharges of retinal ganglion cells deviate systematically from such a relationship. A more realistic relationship can be obtained if the integrate-and-fire mechanism is “leaky”; with this refinement, the integrate-and-fire model captures the essential features of the data. However, the model shows that the distribution of intervals is insensitive to that of the underlying variability. The leakage time constant, threshold, and distribution of the noise are confounded, rendering the model unspecifiable. Another aspect of variability is presented by the variance of responses to repeated discrete stimuli. The variance of response rate increases with the mean response amplitude; the nature of that relationship depends on the duration of the periods in which the response is sampled. These results have defied explanation. But if it is assumed that variability depends on mean rate in the way observed for maintained discharges, the variability of responses to abrupt changes in lighting can be predicted from the observed mean responses. The parameters that provide the best fits for the variability of responses also provide a reasonable fit to the variability of maintained discharges.
INTRODUCTION The axons of the ganglion cells of the vertebrate retina are the sole carriers of visual information from the eyes to the brain. Their message is encoded as a train of action potentials, or neural impulses; it is generally assumed that information is conveyed by the rate of firing impulses. In view of this, it is somewhat enigmatic that the impulse trains are typically quite variable, or “noisy.” The models described
MATHEMATICAL
BIOSCIENCES
112~225-242 (1992)
OElsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010
225 0025-5564/92/$5.00
226
MICHAEL
W. LEVINE
here are concerned with this variability, with an emphasis on how the noise may be introduced. The thesis of this paper is that neural impulses are generated by a mechanism that can be approximated by an integrate-and-fire model [9, 24, 271. In this model, random variability (noise) is added to the deterministic signal conveying the visual information. The combined signal is integrated in time until a threshold is reached; at this point, an impulse is produced, and the integral resets to zero (resting potential). I have been exploring the properties of integrate-and-fire models, using a Monte Carlo digital realization of this scheme (see [14] for an illustration). The static properties of maintained discharges (firing in the absence of any recent change in stimulation) may be characterized by the distribution of the intervals between successive impulses. The properties of this distribution are captured by its moments. Of particular concern are its first moment (the mean inter&, k, whose inverse is the mean firing rate, i.1 and the square root of the second moment about the mean (the standard deuiation, a>. The ratio of the standard deviation to the mean (a / p) is the coefjkient of uariution (CV). This paper will be concerned only with the static properties of the variability, and not with the temporal structure of the noise. There is temporal correlation in the maintained discharges of retinal ganglion cells [e.g., 101, but it is relatively weak in goldfish [131 and even weaker in the cat [6, 23, 241. This temporal correlation is assumed to be associated with the noise [13, 141; the temporal time course of responses to flashes of light is assumed to be a deterministic property of the retinal network. It is likely that some of the deterministic “filtering” mechanisms that act on signals derived from lights also act on the noise to produce autocorrelation in maintained discharges [6, 131. Predictions of integrate-and-fire models will be compared to data from cat ganglion cells, both to test the viability of the model and to see if parameters of the model can be determined from the data. Two kinds of data will be considered: maintained discharges in which the mean firing rate has been set by the presence of steady stimulus lights, and responses to step increments and decrements of contrast.’
‘Contrast generally refers to changes in illumination across a visual display. As used here, it applies to the ratio between the illumination of the center of a receptive field and that of the surround; this may consist of either a centered spot of light superimposed upon a uniform (or dark) field, or a grating of light and dark bars sinusoidally modulated in space. An increment in contrast of a spot of light refers to its onset against the background, and a decrement refers to its offset. An increment of contrast of a grating refers to a 180” phase shift in which a light bar appears on the center; a decrement is a 180” phase shift that places a dark bar on the center.
MODELING
THE VARIABILITY
OF FIRING
RATE
227
METHODS The data have all been reported elsewhere, so a brief description of the methods will suffice. A summary of properties of the cells in these and other relevant studies is presented in Table 1. The studies of maintained discharges [6] were performed at North-
TABLE 1 Comparison
of Properties
in Different
Cells and Studies”
Ganglion cell (retina) Goldfish [191
161
Anesthesia
IsolatedC
Urethane
Maintained discharge’ P (i/s) cv cv vs. 7 Slopeh Log $/a’
26.8 0.89 scattered -0.67 -0.27
Responsej Slope l/4 sk Slope l/2 s’
0.66 0.78
Cat
Cat
30.9 0.69 Figure 1 - 0.55 - 0.03
-
LGNdb
[I61
131
Cat cortex h281
N,O+ halothane
N,O+ halothane
N,O
19.2 0.67
9.4 1.14
low’ 0.6s
- 0.59 - 0.08 0.60 0.71
-0.54 -0.06 0.95 0.91
1.09
“Missing values imply either that the appropriate experiment was not performed or that the result was not reported. bThe LGN, (dorsal lateral geniculate nucleus) is the visual thalamic relay between ganglion cells and cortical cells. ‘Isolated retinas were harvested from euthanized animals. dLGN data were obtained at the same time as the retinal data in [18]; each LGN cell was recorded simultaneously with the ganglion cell that provided its principal input. Reference 4, a reanalysis of data, also supplies some of the data for other columns. e Measures of maintained discharges. ‘Values are not given, but the mean appears to be less than 1 i/s. gValue of CV derived from reported intercept using Eq. (7). hSlope on double logarithmic plot of &,, versus T,. According to Eq. (7), a renewal process should produce a line with a slope of -0.5 [12]. ‘Intercept at T, = 1 s on double logarithmic plot of &,,, versus r,. According to Eq. (7), a renewal process should have an intercept of 0.0 [12]. ‘Measures of variability of responses. kSlope on double logarithmic plot of variance of rate versus rate for 1/4-s samples. ‘Slope on double logarithmic plot of variance of rate versus rate for 1/2-s samples.
228
MICHAEL
W. LEVINE
western University, Evanston, Illinois. Impulses were recorded extracellularly with glass microelectrodes placed in the optic tracts of cats anesthetized with 20-30 mg/(kg-h) urethane. Steady spots slightly larger than the receptive field center were focused on the retina in Maxwellian view. Data were recorded for at least 90 s, commencing after the firing had stabilized. The studies of responses to flashes of light were performed at the University of Sydney, Australia [18]. Impulses were recorded from the ganglion cells of cats anesthetized with a mixture of nitrous oxide (70%) and halothane (l-4%) in oxygen. An intraocular pipette for extracellular recording was placed through a transscleral guide tube. Stimuli were temporally modulated at 1 Hz. Cells were stimulated either by a spot coextensive with the receptive field center, turned on and off at least nine times at each of several luminances, or by an optimal spatial frequency sinusoidal grating, square-wave modulated in counterphase at several contrasts. Responses were taken as the number of impulses produced in a half- or quarter-second period following or preceding any stimulus transition (increase or decrease in luminance at the center of the receptive field). This number was divided by the duration (l/2 or l/4 s) to convert to rate in impulses/second (i/s> for that response time on that trial. Variance of rate was computed from the list of response rates in corresponding times across all the trials. RESULTS
AND
DEPENDENCY
DISCUSSION
OF VARIABILITY
ON MEAN RATE
Earlier versions of the digital integrate-and-fire model used a simple cumulation of the sum of the rate-setting signal (a constant referred to as bias) and the random noise [14]. Mean rate was shifted by selecting different values of the bias; it was observed that the relationship of CV to mean rate was an inverse square root; CV=
ar-‘/*,
(1)
where a is a constant. This relationship could have been derived by noting that the simple integrate-and-fire is equivalent to a random walk with superimposed drift [3, 7, 8, 221. The diffusion approximation to that model provides equations for the moments of the predicted distribution that predict the inverse square root relationship assuming that mean rate is shifted by altering the “drift rate” parameter (see [151X As a first approximation, the maintained discharges comply with the inverse square root relationship. The logarithm of CV was plotted versus the logarithm of the mean rate (for each cell), and the best-fit (principal-component) line was calculated. The mean of the slopes of
MODELING
THE VARIABILITY
-
OF FIRING
RATE
229
-025
2 B -.
-050
-100'
’
00
05
’
10
’
15
’ 20
’
25
Log (rate)
FIG. 1. Relationship between CV and mean rate (impulses/s), on double logarithmic axes. Solid stars represent averaged data from 10 retinal ganglion cells in the cat [6], averages taken every 0.1 log unit on the log rate axis. Circles are derived from simulated firing generated by a leaky integrate-and-fire model (10 ms time constant). The smooth curve is the average relationship from models in which the CV versus rate relationship was used to predict variability of rate [17], fit to data from 13 cat ganglion cells [18].
these lines is -0.51 (N = 10) 1141. While an inverse square root relationship is the best power law relationship, it is clearly not the optimal fit. The data from individual cells show a distinct downward curvature on these plots [14], which can also be seen on the plot of the average log CV versus log mean rate, shown as the stars in Figure 1. In this graph, all the available cat maintained discharge data were averaged in increments of 0.1 log unit (of rate). The first attempt to explain this discrepancy was based on the observation that the effect of threshold level is opposite to that of bias [14]. If the mean rate is shifted by an adjustment of threshold, CV is directly proportional to the square root of the mean rate. This is clearly contrary to the observations, but it is possible to imagine an adjustment of threshold that would affect the shape of the curve. In particular, an exponential refractory period, in which threshold becomes temporarily higher immediately following each impulse and then decays to its steady value, produces some curvature in the appropriate direction [14]. However, this modification generally could not produce enough curvature over an extended domain to explain the observed relationship.
MICHAEL
230
W. LEVINE
No such modification is required if an inherent absurdity of the simple model is removed. In the model discussed so far, the integration has perfect memory. Thus, if the input signal (bias plus noise) were suddenly removed (set to zero> just before threshold was attained, the integral would hold its near-threshold value indefinitely; when the input was restored, the integrator would have only to cumulate the short distance from that previous level to threshold to produce an impulse. Real neurons do not have this property: in the absence of an input, the membrane potential decays to resting potential (defined as zero in the model). That is, neurons may be said to be leaky [9]. Leakage is not the same as drift, because the absolute rate of change (mV/s> is greater when the membrane potential is farther from resting potential. Its effect (like that of the “refractory” threshold) is thus different at high rates, where potential rapidly rises to threshold so inputs cumulate, and low rates, where the potential decays in the intervals between impulses so threshold is reached by individual large input values.2 Leakage is implemented in the model as exponential decay, with the time constant as an additional parameter. The value of the integral at each millisecond is its previous value plus the input (bias plus noise); if that value does not achieve threshold, it is multiplied by e-‘I’, the decay that would occur in the 1 ms until the succeeding input is added CT is the time constant). The effect of leakage may be seen in Figure 2, a comparison of models with various time constants. Time constants in a range appropriate for the physiology of neurons produce curvature like that observed in the ganglion cell data. This is demonstrated by the open circles in Figure 1, which were derived from a model with a 10 ms time constant; the model fits the mean data almost perfectly. Although the curves in Figure 2 demonstrate that the leaky integrate-and-fire model reproduces the observed effect, it unfortunately is not true that we may now use those data to specify the leakage time constant of the ganglion cells. There is a confounding of time constant and threshold that allows models with the same time constant to display different curvatures. This may be seen in Figure 3. In this figure, all the models have 5 ms time constants; they differ in threshold. The general form of the leaky integrate-and-fire model is what is required to fit the data, but the models are not unique. INTERVAL
DISTRIBUTIONS
One advantage of the digital realization of the integrate-and-fire model is that it is relatively simple to tailor the distribution of the noise
‘Since the inputs do not cumulate at low rates, the firing approaches process, and CV + 1. This effect may be observed in Figures l-3.
a Poisson
MODELING
THE VARIABILITY
OF FIRING
231
RATE
_....
\
\
I
I
I
10
05
vu
Log
\
\
j
I
15
\
20
25
(rate)
FIG. 2. Comparison of relationships between CV and mean rate for integrateand-fire models with various leakage time constants. Solid squares show two models with no leakage (T =m); the two symbol sizes represent two different thresholds. Other time constants: asterisks = 50 ms; diamonds = 20 ms; circles = 10 ms; stars = 5 ms; triangles = 2 ms. Uniform distribution noise for all models. Dashed lines have a slope of -0.5.
025 5 000
_
0 o.**.*~.o...o~
._........
8; 5
-025
._......
* 0
-
_._ __....._
Oo * * *
B -J
0 -050
*
0 0
-075
-100 00
,
I
I
I
05
10
15
20
25
Log (rate)
FIG. 3. Comparison of relationships between CV and mean rate for integrateand-fire models with various thresholds but the same time constant (5 msl. Threshold values: 0.5 (circles), 0.2 (stars), and 5.0 (diamonds). Gaussian noise for all models.
232
MICHAEL
W. LEVINE
that is added to the bias. In most analytical treatments, it would be difficult to evaluate the effects of different distributions that might be assumed for various theoretical reasons. One would expect the form of the noise distribution to have its principal effect on the shape of the interval distribution; indeed, some workers have used the shape of the interval distribution to attempt to establish the source of the noise (e.g., [2]). To test the validity of that idea, models were generated in which the noise was Gaussian (zero skew, normokurtic), uniformly distributed (zero skew, platykurtic), or first-order gamma (positive skew, leptokurtic); the interval distributions generated with these noise distributions were indistinguishable when bias (and/or leakage) was adjusted to give approximately equal mean rate and CV [15]. The examples in Figure 4 show models that simulate the firing normally observed in cat (top row) and firing with a relatively low CV (bottom row). The distributions in Figure 4 are fit with two theoretical curves: the log normal distribution and the diffusion approximation to the random walk with drift. These distributions generally fit better than others that have been used for interval distributions E1.51,including the hyperbolic normal distribution [161, the hyperbolic gamma distribution [25], and the gamma distribution [l, 10,231. The random walk approximation typically fits the simulated data slightly better than the log normal distribution, which is perhaps to be expected because the model (in the absence of leakage) is equivalent to the random walk with drift model. However, the random walk approximation also works well when the integrateand-fire is made leaky. Note that the fits shown for each row of Figure 4 are identical; instead of the best fits to each, the mean best fit is used for all three distributions with a similar CV. Figure 5 shows that the interval distributions generated by the model are appropriate for data. These interval distributions from cat ganglion cells are strikingly similar to the simulations in Figure 4. The log normal distributions and random walk approximations shown on these distributions are identical to those used in Figure 4. In fits to data, the log normal distribution often provides a better fit than the random walk approximation; it also is more forgiving of irregularities such as a possible weak second mode in the interval distribution. Given the ambiguity about the true nature of the distribution, it is probably wiser to use the log normal distribution, for it does not imply particular parameters of the spike-generating process [ 151. While it appears that the interval distributions do not reflect the “signature” of the noise distribution used to generate them, it remains possible that some trace of that distribution can be found in higher moments of the interval distributions. The third moment about the
MQDELING THE VARIABILITY OF FIRfNG RATE
233
234
MICHAEL
Interval
W. LEVINE
(ms)
FIG. 5. Interval distributions from two cat ganglion cells. The smooth curves are identical to those in the corresponding rows of Figure 4. Top: Off-X cell, CV = 0.76 at 21.6 i/s. Bottom: On-X cell, CV = 0.37 at 61.8 i/s. Both recorded from optic tract
161.
mean
is the skew, S:
s=/=(T-?)‘p(+k
(2)
--o
Since P(T) = 0 for T =G0, we may change with some rearranging and substitution, s=
the lower integration
llP3PWd~ 03
1 --
--
(y”
limit, and
3 CV’
Substitution of each of the putative distributions as ~(7) in Eq. (3) gives relationships that can be solved for skew in terms of CV [26]. Specifically, for the random walk approximation (as given in [3]), S random walk= For the log normal
distribution S lognormal
3 cv.
(4)
(as given in [15]), =3cv+cv”.
(5)
MODELING
Finally,
THE VARIABILITY
for the hyperbolic
OF FIRING
gamma
shyperbolic
distribution
RATE
235
(suggested
4cv I- =
in [25]),
(6)
~
l-cv*’
These relationships could allow an additional test of the suggested distributions. Notice that the hyperbolic gamma distribution, which in any case provides a poorer fit than the other two, incorporates the absurdity of S + cc as CV + 1, and S < 0 for CV > 1. Skew is plotted as a function of CV in Figure 6. Various models, using different noise distributions, some with leakage and some without, are shown; for each model, CV (and rate) was adjusted by changing bias. These models produce skews that cluster near and below the 3 CV line predicted for the random walk approximation. There is no obvious discrimination according to leakage or threshold. Data from ganglion cells also fall in a similar band on this plot. The three cells from which data were obtained at the widest range of firing rates (and hence of CV) are shown by solid symbols in Figure 6. Data from other ganglion cells [18] fall in a similar band (not shown on the
0.2
0.4
0.6
0.0
1.0
1.2
FIG. 6. Relationship between skew and CV derived from various integrate-and-fire simulations and data from cat. Simulations: asterisks = Gaussian noise distribution, no leakage; stars = first-order gamma noise, no leakage. Other symbols for uniform noise distributions (square = no leakage; circle = 5 ms time constant, triangles = 2 ms time constant, with two different threshold levels represented by erect or inverted triangles). Solid symbols are data from three cat ganglion cells [6]; solid circles = OffX, solid squares = Off-X, solid diamonds = On-Y.
236
MICHAEL
W. LEVINE
plot). For some cells, the data were noticeably lower than the 3 CV line; these cells display a slight tendency toward bimodal interval distributions. A bimodal distribution can have a relatively larger standard deviation and smaller skew. The models all produce skews comparable to those of the data; however, there is a tendency for the Gaussian noise distribution to produce somewhat higher skews, particularly at high CV. (The Gaussian simulations can produce even higher skews than the ones shown in Figure 6.) Higher skew could be due to the extended lower tail (negative values) of the Gaussian distribution, which both the gamma and uniform distributions lack. This would appear to disqualify Gaussian noise, except that the differences are not significant for the number of simulations and data so far analyzed. If the difference is real, it can easily be ameliorated by postulating a weakly rectifying nonlinearity that attenuates the lower tail of the Gaussian noise. While skew seems not to discriminate among models, it may seem to favor the random walk approximation for fitting interval distributions. However, the relatively low skew (at high CV) is an artifact of the limited data available from a neuron. The theoretical curves were derived by integration to infinity; real data do not include those exquisitely rare very long intervals that are nonetheless inherent in the theoretical curve. Indeed, if an interval in excess of several seconds did occur, the experimenter would assume that some mishap had befallen and discard the run. Because of the cubing operation [see Eq. (211, these rare long intervals have a large effect on skew. Numerical integration of the moments of the log normal distribution shows that for CV near 1 (the fit shown in Figure 2b in [15]), skew is less than 75% of its expected value (for the measured CV) when the area is more than 0.999 (first 325 ms). Thus, the predicted skew for a real interval distribution that is a perfect log normal would be the same as that for the theoretical random walk approximation. This effect is most pronounced for high CV, where high skew is expected, and renders the comparisons in Figure 6 useless. Obviously, examination of higher moments would be futile, for they will be even more disrupted. VARLABILITY
OF RESPONSES
Another way of looking at variability was introduced by workers examining responses of cat cortical cells [5, 28, 291. Cortical cells have little or no maintained discharges, but the responses to repeated stimuli are quite variable. Despite considerable scatter in plots of log variance of rate versus log rate, Tolhurst and his coworkers [5, 28, 291 were able to show that the variance of rate of responses was directly proportional to the mean response amplitude (rate). This was later shown to also be
MODELING
THE VARIABILITY
OF FIRING
RATE
237
true for lateral geniculate neurons (the thalarnic cells that precede cortical cells in the signal pathway from retina to brain) [4]. This result is inconsistent with the nearly inverse square root relationship between CV and mean rate described for cat ganglion cells 114, 231, embodied in Eq. (1). If the firing is a renewal process, as it very nearly is in cat ganglion cells [6], then
where T, is the duration of the samples in which rate is measured [12]. Substitution of Eq. (1) into Eq. (7) shows that variance of rate should be independent of mean rate for any given sample duration [17, 191. It was therefore suggested that the dependence of variance of rate upon mean rate in cortical cells represents a difference between ganglion cells and cortical cells. To test whether there really is a difference between cortical cells and ganglion cells, Levine et al. [191 measured the variance of rate of goldfish ganglion cells in responses to repeated stimuli. They found that variance of rate clearly increased as a function of mean rate, contrary to the expectation from maintained discharges (of cat ganglion cells). Moreover, since ganglion cells support a moderate mean maintained discharge rate, variance of firing was also measured during the phases of the responses that represent decreases in firing rate; variance of rate was a function of firing rate during the response, and not of the strength of the stimulus that elicited that response. There were also some differences from the findings from cortical cells: the variance of rate increased less than in direct proportion to the rate, and the exact nonlinear relationship between variance of rate and rate depended on the length of the periods over which the responses were sampled. In particular, a plot of the logarithm of variance of rate versus the logarithm of mean rate showed that the data were reasonably well fit by a straight line with a fractional slope. The longer the period after each stimulus onset or offset in which responses were measured, the steeper the slope of the best-fit line (see Table 1). This was particularly puzzling: Why would responses measured for longer periods (which thereby include more of the steady maintained-like “plateau”) show more nonlinear interaction with noise than those measured in brief periods (the “purer” responses)? If the variance of rate of maintained discharges is independent of rate [Eq. (7)], the longer durations should show less dependence on mean rate, not more. Since the measurements of variance of rate of ganglion cells were
238
MICHAEL W. LEVINE
made using a different experimental animal, the results did not directly address the question of whether the apparent contradiction between the outcome of the cortical experiments on repeated stimuli and the ganglion cell experiments on maintained discharges was a difference between retinal versus cortical cells or a difference between responses to abrupt stimulation versus maintained discharges. To settle this issue, similar experiments on variability of responses were performed on cat ganglion cells; the results were nearly identical to those from the fish [18] (see Table 1). The only deviation noted was that about one-third of the cells tested showed a distinctly curved relationship on the log variance versus log rate plots. The species did not account for the contradiction; the difference appears to be between responses to sudden transitions in lighting versus maintained discharges in temporally steady conditions. The apparent contradiction between the experiments on variance of rate and those on maintained discharges raised the specter of two kinds of variability that are handled differently in the retina [19]. The noise that affects maintained discharges is added to the rate-setting signal (bias); the noise responsible for the variability of responses to abrupt transitions seems to have a nonlinear relationship with the rate-setting mechanism. How might the retinal network add noise to the steady signal setting maintained rate while multiplying the noise by the transient signals effecting responses ? And, as asked above, why is the multiplication more evident for longer sample durations, which should be more like the maintained discharges? Recent work has shown that these results can be reconciled [17]. As the nonlinearity of the double logarithmic plot derived from steady illumination (Figure 1) illustrates, CV is not actually proportional to the inverse square root of mean rate. Moreover, the firing rate during a response to an abrupt change in illumination of the retina extends through a range of values that may be quite different from the mean rate during the period in which the response is measured. Typically, there is a peak (increase or decrease) immediately upon stimulation, followed by a gradual relaxation to a plateau nearer the maintained firing rate. The longer the duration of the measurement period, the more different the extreme peak will be from the overall mean. The stronger the stimulus, the more extreme the difference will be and the more sharply the peak will decline to the plateau. Levine and Zimmerman [17] proposed a model that, by taking these facts into account, reconciles the two kinds of data. The premise of our model is that the same additive noise affects responses to sudden changes as produces the variability of maintained discharges. The variability is represented by CV, which was assumed to
MODELING
THE VARIABILITY
OF FIRING
RATE
change as a function of rate (F) in the way shown in Figure specifically, a parabolic relationship was assumed:
239
1;
The coefficient of variation is computed from the responses in the response histogram, which is measured in short bins (synchronized to the stimulus and averaged over all the repetitions of the stimulus). This histogram captures the time course of the response, tracing its evolution from peak to plateau. One could not measure CV in such short periods (during a single cycle of which only one or two impulses might occur), but from the mean rate and Eq. (8) one can infer the CV that would be expected. From CV, one can calculate the variance of count expected in that bin; assuming independence of the bins that comprise a response, the variances of count simply sum, and variance of rate can be found by dividing by duration squared. Thus, the variance of rate derived from the mean responses in short time periods can be compared to the observed variance of rate. The model was fit to the data from ganglion cells [18, 191, with the coefficients of the parabolic fit of log CV versus log rate as the only free parameters.3 The model was able to fit the data at least as well as linear regression (despite having fewer free parameters when fitting data from goldfish), even reproducing some of the “scatter” evident in the plots of log variance of rate versus log rate. In those cases in which cat ganglion cells produced curved relationships, the model reproduced the curvature [17]. The mean parabolic relationship from the cat models is plotted as a solid line in Figure 1. It is of the right general form of fit the observed maintained discharges, which are well fit by the leaky integrate-and-fire model. The vertical error is not statistically significant and is necessarily linked to the actual CV of the sample (see footnote 3). Remember that the maintained discharge data and the data upon which the response
3Equation (8) has three coefficients, but only two are used as parameters of the model. A is found from X, and .r2 given the mean rate and CV of the observed maintained discharge of the cell. The third model parameter is a multiplier used in converting CV to variance of rate; this fractional factor (very nearly unity) accounts for the slight temporal structure of the firing, which is not exactly a renewal process 16, 12, 131. The value of the third parameter is consistent with the values in rows 8 and 9 of Table 1. Since multiplication of CV is equivalent to adding a constant to log CV, this is equivalent to allowing adjustment of A. However, the values of A used for the smooth curve in Figure 1 remain set to reproduce the observed CV of each cell.
240
MICHAEL
W. LEVINE
variability model were based were drawn from experiments, in laboratories hemispheres apart, thetics.
two different sets of using different anes-
CONCLUSION Our basic model suggests that the variable discharges of a ganglion cell are generated by adding noise to the deterministic signal representing visual information, and using this combined signal as the input to a leaky integrate-and-fire mechanism. This model accounts well for the variability of maintained discharges, both in terms of the distribution between successive neural impulses and in the way in which variability changes as mean rate shifts in response to steady stimulation. It also can reproduce the higher order statistics and temporal dependencies in the firing [6, 13, 141. It is also consistent with the variability of firing in response to transient stimulation. data, it must be Although the model accounts for “well-behaved” noted that interval distributions from ganglion cells are not always smoothly unimodal. A minority of ganglion cells produce multimodal interval distributions, and these cannot always be dismissed as pathological (as some workers do). The models considered here make no attempt to explain multimodal distributions, although a modification in which the threshold changes after each impulse can produce such an effect [14]. Alternatively, multimodal distributions could be due to a “multiplexing” among different noise or signal sources [16] or to a second mechanism by which impulses may be produced according to the derivative of the membrane potential [ll, 20, 211. These complications have not been considered here. Finally, it is important to remember that even if an integrate-and-fire model is satisfactory, it is not unique. The trade-offs among parameters allow the same data to be fit with a large number of models; intracellular measurements will be required to determine the actual parameters of ganglion cells. I am gratefil to the many collaborators over the years who have offered advice and suggestions as this work progressed; major contributions were made by Violeta Cam’&-Carire, Jeremy Shefner, and especially Roger Zimmerman. I am particularly grateful to those collaborators who have allowed me to share their facilities and harness their expertise to collect the data reported here, speci@ally Laura Frishman and Christina Enroth-Cugell at Northwestern University and Brian Cleland at the University of Sydney. I also wish to express my appreciation to Grace Yang and Charles Smith, organizers of this conference, for inviting me to participate.
MODELING
THE VARIABILITY
OF FIRING
RATE
241
REFERENCES
4
10 11 12 13
H. B. Barlow and W. R. Levick, Changes in the maintained discharge with adaptation level in the cat retina J. Physiol. 202:699-718 (1969). H. B. Barlow, W. R. Levick, and M. Yoon, Responses to single quanta of light in retinal ganglion cells of the cat, I/is&n Res. ll(Supp1. 3):87-101 (1971). D. Berger, K. Pribram, H. Wild, and C. Bridges, An analysis of neural spike-train distributions: determinants of the response of visual cortex neurons to changes in orientation and spatial frequency, Exp. Brain. Res. 80:129-134 (1990). V. Carrion-Carire, B. G. Cleland, A. W. Freeman, M. W. Levine, and R. P. Zimmerman, Variability of responses in retinal and geniculate cells, Int,est. Ophthalmol. Vis. Sci. 29(Suppl):295 (1988) (Abstract). A. F. Dean, The variability of discharge of simple cells in the cat striate cortex, Exp. Bruin Res. 44:437-440 (1981). L. J. Frishman and M. W. Levine, Statistics of the maintained discharge of cat retinal ganglion cells, J. Physiol. 339:475-494 (1983). G. L. Gerstein and B. Mandelbrot, Random walk models for the spike activity of a single neuron, Biophys. J. 4:41-68 (1964). P. I. M. Johannesma, Diffusion models for the stochastic activity of neurons, in Neural Networks, E. R. Caianiello, Ed., Springer-Verlag, Berlin, 1968, pp. 116-144. B. W. Knight, Dynamics of encoding in a population of neurons, J. Gen. Physiol. 591734-766 (1972). S. W. Kuffler, R. FitzHugh,
and H. B. Barlow, Maintained activity in the cat’s retina in light and darkness, J. Gen. Physiol. 40:683-702 (1957). M. J. M. Lankheet, J. Molenaar, and W. A. van de Grind, The spike generating mechanism of cat retinal ganglion cells, Vision Res. 29:.505-517 (1989). M. W. Levine, Firing rates of a retinal neuron are not predictable from interspike interval statistics, Biophys. /. 30:9-26 (1980). M. W. Levine, Retinal processing of intrinsic and extrinsic noise, J. Neurophysiol. 48:992-1010
(1982).
M. W. Levine, Variability in the maintained discharges of retinal ganglion cells, J. Opt. Sot. Am. A 4:2308-2320 (1987). I5 M. W. Levine, The distribution of the intervals between neural impulses in the maintained discharges of retinal ganglion cells, Biol. Cybem. 65:459-467 (1991). 16 M. W. Levine and J. M. Shefner, A model for the variability of interspike intervals during sustained firing of a retinal neuron, Biophys. J. 19:241-252 14
(1977).
17 18 19 20
M. W. Levine and R. P. Zimmerman, A model for the variability of maintained discharges and responses to flashes of light, Biol. Cybem. 65~469-477 (1991). M. W. Levine, B. G. Cleland, and R. P. Zimmerman, Variability of responses of cat retinal ganglion cells, Visual Neurosci. 8:277-279 (1992). M. W. Levine, R. P. Zimmerman, and V. Carrion-Carire, Variability in responses of retinal ganglion cells, /. Opt. Sot. Am. A 5:593-597 (1988). R. Llinas and H. Jahnsen, Electrophysiology of mammalian thalamic neurones in vitro, Nature 297:406-408
21
(1982).
F.-S. Lo and S. M. Sherman, In vivo recording of postsynaptic potentials and low threshold spikes in W cells of the cat’s lateral geniculate nucleus, Exp. Bruin Res. 81:438-442 (1990).
242 22 23 24 25 26 27 28
29
MICHAEL
W. LEVINE
J. Pernier and P. Gerin, Temporal pattern analysis of spontaneous unit activity in the neocortex, Biol. Cybem. 18:123-136 (19751. J. G. Robson and J. B. Troy, The nature of the noise in cat retinal ganglion cells, J. Opt. Sot. Am. A 4:2301-2307 (1987). R. W. Rodieck, Maintained activity of cat retinal ganglion cells, J. Neurophysiol. 30:1043-1071 (1967). C. E. Smith, A comment on a retinal neuron model, Biophys. /. 25:385-386 (1979). C. E. Smith, personal communication (1991). R. B. Stein, A theoretical analysis of neuronal variability, Biophys. J. 5:173-194 (1965). D. J. Tolhurst, J. A. Movshon, and A. F. Dean, The statistical reliability of signals in single neurons in cat and monkey visual cortex, Ksion Res. 23:775-785 (1983). D. J. Tolhurst, J. A. Movshon, and I. D. Thompson, The dependence of response amplitude and variance of cat visual cortical neurones on stimulus contrast, Exp. Bruin Res. 41:414-419 (1981).
W. LEVINE
Uniuersityof Illinois at Chicago, Department of Psychology and Committee on Neuroscience,
Chicago, Illinois 60680 Received 4 October 1991; reuised 28 July 1992
ABSTRACT Impulse trains simulating the maintained discharges of retinal ganglion cells were generated by digital realizations of the integrate-and-fire model. If the mean rate were set by a “bias” level added to “noise,” the variability of firing would be related to the mean firing rate as an inverse square root law; the maintained discharges of retinal ganglion cells deviate systematically from such a relationship. A more realistic relationship can be obtained if the integrate-and-fire mechanism is “leaky”; with this refinement, the integrate-and-fire model captures the essential features of the data. However, the model shows that the distribution of intervals is insensitive to that of the underlying variability. The leakage time constant, threshold, and distribution of the noise are confounded, rendering the model unspecifiable. Another aspect of variability is presented by the variance of responses to repeated discrete stimuli. The variance of response rate increases with the mean response amplitude; the nature of that relationship depends on the duration of the periods in which the response is sampled. These results have defied explanation. But if it is assumed that variability depends on mean rate in the way observed for maintained discharges, the variability of responses to abrupt changes in lighting can be predicted from the observed mean responses. The parameters that provide the best fits for the variability of responses also provide a reasonable fit to the variability of maintained discharges.
INTRODUCTION The axons of the ganglion cells of the vertebrate retina are the sole carriers of visual information from the eyes to the brain. Their message is encoded as a train of action potentials, or neural impulses; it is generally assumed that information is conveyed by the rate of firing impulses. In view of this, it is somewhat enigmatic that the impulse trains are typically quite variable, or “noisy.” The models described
MATHEMATICAL
BIOSCIENCES
112~225-242 (1992)
OElsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010
225 0025-5564/92/$5.00
226
MICHAEL
W. LEVINE
here are concerned with this variability, with an emphasis on how the noise may be introduced. The thesis of this paper is that neural impulses are generated by a mechanism that can be approximated by an integrate-and-fire model [9, 24, 271. In this model, random variability (noise) is added to the deterministic signal conveying the visual information. The combined signal is integrated in time until a threshold is reached; at this point, an impulse is produced, and the integral resets to zero (resting potential). I have been exploring the properties of integrate-and-fire models, using a Monte Carlo digital realization of this scheme (see [14] for an illustration). The static properties of maintained discharges (firing in the absence of any recent change in stimulation) may be characterized by the distribution of the intervals between successive impulses. The properties of this distribution are captured by its moments. Of particular concern are its first moment (the mean inter&, k, whose inverse is the mean firing rate, i.1 and the square root of the second moment about the mean (the standard deuiation, a>. The ratio of the standard deviation to the mean (a / p) is the coefjkient of uariution (CV). This paper will be concerned only with the static properties of the variability, and not with the temporal structure of the noise. There is temporal correlation in the maintained discharges of retinal ganglion cells [e.g., 101, but it is relatively weak in goldfish [131 and even weaker in the cat [6, 23, 241. This temporal correlation is assumed to be associated with the noise [13, 141; the temporal time course of responses to flashes of light is assumed to be a deterministic property of the retinal network. It is likely that some of the deterministic “filtering” mechanisms that act on signals derived from lights also act on the noise to produce autocorrelation in maintained discharges [6, 131. Predictions of integrate-and-fire models will be compared to data from cat ganglion cells, both to test the viability of the model and to see if parameters of the model can be determined from the data. Two kinds of data will be considered: maintained discharges in which the mean firing rate has been set by the presence of steady stimulus lights, and responses to step increments and decrements of contrast.’
‘Contrast generally refers to changes in illumination across a visual display. As used here, it applies to the ratio between the illumination of the center of a receptive field and that of the surround; this may consist of either a centered spot of light superimposed upon a uniform (or dark) field, or a grating of light and dark bars sinusoidally modulated in space. An increment in contrast of a spot of light refers to its onset against the background, and a decrement refers to its offset. An increment of contrast of a grating refers to a 180” phase shift in which a light bar appears on the center; a decrement is a 180” phase shift that places a dark bar on the center.
MODELING
THE VARIABILITY
OF FIRING
RATE
227
METHODS The data have all been reported elsewhere, so a brief description of the methods will suffice. A summary of properties of the cells in these and other relevant studies is presented in Table 1. The studies of maintained discharges [6] were performed at North-
TABLE 1 Comparison
of Properties
in Different
Cells and Studies”
Ganglion cell (retina) Goldfish [191
161
Anesthesia
IsolatedC
Urethane
Maintained discharge’ P (i/s) cv cv vs. 7 Slopeh Log $/a’
26.8 0.89 scattered -0.67 -0.27
Responsej Slope l/4 sk Slope l/2 s’
0.66 0.78
Cat
Cat
30.9 0.69 Figure 1 - 0.55 - 0.03
-
LGNdb
[I61
131
Cat cortex h281
N,O+ halothane
N,O+ halothane
N,O
19.2 0.67
9.4 1.14
low’ 0.6s
- 0.59 - 0.08 0.60 0.71
-0.54 -0.06 0.95 0.91
1.09
“Missing values imply either that the appropriate experiment was not performed or that the result was not reported. bThe LGN, (dorsal lateral geniculate nucleus) is the visual thalamic relay between ganglion cells and cortical cells. ‘Isolated retinas were harvested from euthanized animals. dLGN data were obtained at the same time as the retinal data in [18]; each LGN cell was recorded simultaneously with the ganglion cell that provided its principal input. Reference 4, a reanalysis of data, also supplies some of the data for other columns. e Measures of maintained discharges. ‘Values are not given, but the mean appears to be less than 1 i/s. gValue of CV derived from reported intercept using Eq. (7). hSlope on double logarithmic plot of &,, versus T,. According to Eq. (7), a renewal process should produce a line with a slope of -0.5 [12]. ‘Intercept at T, = 1 s on double logarithmic plot of &,,, versus r,. According to Eq. (7), a renewal process should have an intercept of 0.0 [12]. ‘Measures of variability of responses. kSlope on double logarithmic plot of variance of rate versus rate for 1/4-s samples. ‘Slope on double logarithmic plot of variance of rate versus rate for 1/2-s samples.
228
MICHAEL
W. LEVINE
western University, Evanston, Illinois. Impulses were recorded extracellularly with glass microelectrodes placed in the optic tracts of cats anesthetized with 20-30 mg/(kg-h) urethane. Steady spots slightly larger than the receptive field center were focused on the retina in Maxwellian view. Data were recorded for at least 90 s, commencing after the firing had stabilized. The studies of responses to flashes of light were performed at the University of Sydney, Australia [18]. Impulses were recorded from the ganglion cells of cats anesthetized with a mixture of nitrous oxide (70%) and halothane (l-4%) in oxygen. An intraocular pipette for extracellular recording was placed through a transscleral guide tube. Stimuli were temporally modulated at 1 Hz. Cells were stimulated either by a spot coextensive with the receptive field center, turned on and off at least nine times at each of several luminances, or by an optimal spatial frequency sinusoidal grating, square-wave modulated in counterphase at several contrasts. Responses were taken as the number of impulses produced in a half- or quarter-second period following or preceding any stimulus transition (increase or decrease in luminance at the center of the receptive field). This number was divided by the duration (l/2 or l/4 s) to convert to rate in impulses/second (i/s> for that response time on that trial. Variance of rate was computed from the list of response rates in corresponding times across all the trials. RESULTS
AND
DEPENDENCY
DISCUSSION
OF VARIABILITY
ON MEAN RATE
Earlier versions of the digital integrate-and-fire model used a simple cumulation of the sum of the rate-setting signal (a constant referred to as bias) and the random noise [14]. Mean rate was shifted by selecting different values of the bias; it was observed that the relationship of CV to mean rate was an inverse square root; CV=
ar-‘/*,
(1)
where a is a constant. This relationship could have been derived by noting that the simple integrate-and-fire is equivalent to a random walk with superimposed drift [3, 7, 8, 221. The diffusion approximation to that model provides equations for the moments of the predicted distribution that predict the inverse square root relationship assuming that mean rate is shifted by altering the “drift rate” parameter (see [151X As a first approximation, the maintained discharges comply with the inverse square root relationship. The logarithm of CV was plotted versus the logarithm of the mean rate (for each cell), and the best-fit (principal-component) line was calculated. The mean of the slopes of
MODELING
THE VARIABILITY
-
OF FIRING
RATE
229
-025
2 B -.
-050
-100'
’
00
05
’
10
’
15
’ 20
’
25
Log (rate)
FIG. 1. Relationship between CV and mean rate (impulses/s), on double logarithmic axes. Solid stars represent averaged data from 10 retinal ganglion cells in the cat [6], averages taken every 0.1 log unit on the log rate axis. Circles are derived from simulated firing generated by a leaky integrate-and-fire model (10 ms time constant). The smooth curve is the average relationship from models in which the CV versus rate relationship was used to predict variability of rate [17], fit to data from 13 cat ganglion cells [18].
these lines is -0.51 (N = 10) 1141. While an inverse square root relationship is the best power law relationship, it is clearly not the optimal fit. The data from individual cells show a distinct downward curvature on these plots [14], which can also be seen on the plot of the average log CV versus log mean rate, shown as the stars in Figure 1. In this graph, all the available cat maintained discharge data were averaged in increments of 0.1 log unit (of rate). The first attempt to explain this discrepancy was based on the observation that the effect of threshold level is opposite to that of bias [14]. If the mean rate is shifted by an adjustment of threshold, CV is directly proportional to the square root of the mean rate. This is clearly contrary to the observations, but it is possible to imagine an adjustment of threshold that would affect the shape of the curve. In particular, an exponential refractory period, in which threshold becomes temporarily higher immediately following each impulse and then decays to its steady value, produces some curvature in the appropriate direction [14]. However, this modification generally could not produce enough curvature over an extended domain to explain the observed relationship.
MICHAEL
230
W. LEVINE
No such modification is required if an inherent absurdity of the simple model is removed. In the model discussed so far, the integration has perfect memory. Thus, if the input signal (bias plus noise) were suddenly removed (set to zero> just before threshold was attained, the integral would hold its near-threshold value indefinitely; when the input was restored, the integrator would have only to cumulate the short distance from that previous level to threshold to produce an impulse. Real neurons do not have this property: in the absence of an input, the membrane potential decays to resting potential (defined as zero in the model). That is, neurons may be said to be leaky [9]. Leakage is not the same as drift, because the absolute rate of change (mV/s> is greater when the membrane potential is farther from resting potential. Its effect (like that of the “refractory” threshold) is thus different at high rates, where potential rapidly rises to threshold so inputs cumulate, and low rates, where the potential decays in the intervals between impulses so threshold is reached by individual large input values.2 Leakage is implemented in the model as exponential decay, with the time constant as an additional parameter. The value of the integral at each millisecond is its previous value plus the input (bias plus noise); if that value does not achieve threshold, it is multiplied by e-‘I’, the decay that would occur in the 1 ms until the succeeding input is added CT is the time constant). The effect of leakage may be seen in Figure 2, a comparison of models with various time constants. Time constants in a range appropriate for the physiology of neurons produce curvature like that observed in the ganglion cell data. This is demonstrated by the open circles in Figure 1, which were derived from a model with a 10 ms time constant; the model fits the mean data almost perfectly. Although the curves in Figure 2 demonstrate that the leaky integrate-and-fire model reproduces the observed effect, it unfortunately is not true that we may now use those data to specify the leakage time constant of the ganglion cells. There is a confounding of time constant and threshold that allows models with the same time constant to display different curvatures. This may be seen in Figure 3. In this figure, all the models have 5 ms time constants; they differ in threshold. The general form of the leaky integrate-and-fire model is what is required to fit the data, but the models are not unique. INTERVAL
DISTRIBUTIONS
One advantage of the digital realization of the integrate-and-fire model is that it is relatively simple to tailor the distribution of the noise
‘Since the inputs do not cumulate at low rates, the firing approaches process, and CV + 1. This effect may be observed in Figures l-3.
a Poisson
MODELING
THE VARIABILITY
OF FIRING
231
RATE
_....
\
\
I
I
I
10
05
vu
Log
\
\
j
I
15
\
20
25
(rate)
FIG. 2. Comparison of relationships between CV and mean rate for integrateand-fire models with various leakage time constants. Solid squares show two models with no leakage (T =m); the two symbol sizes represent two different thresholds. Other time constants: asterisks = 50 ms; diamonds = 20 ms; circles = 10 ms; stars = 5 ms; triangles = 2 ms. Uniform distribution noise for all models. Dashed lines have a slope of -0.5.
025 5 000
_
0 o.**.*~.o...o~
._........
8; 5
-025
._......
* 0
-
_._ __....._
Oo * * *
B -J
0 -050
*
0 0
-075
-100 00
,
I
I
I
05
10
15
20
25
Log (rate)
FIG. 3. Comparison of relationships between CV and mean rate for integrateand-fire models with various thresholds but the same time constant (5 msl. Threshold values: 0.5 (circles), 0.2 (stars), and 5.0 (diamonds). Gaussian noise for all models.
232
MICHAEL
W. LEVINE
that is added to the bias. In most analytical treatments, it would be difficult to evaluate the effects of different distributions that might be assumed for various theoretical reasons. One would expect the form of the noise distribution to have its principal effect on the shape of the interval distribution; indeed, some workers have used the shape of the interval distribution to attempt to establish the source of the noise (e.g., [2]). To test the validity of that idea, models were generated in which the noise was Gaussian (zero skew, normokurtic), uniformly distributed (zero skew, platykurtic), or first-order gamma (positive skew, leptokurtic); the interval distributions generated with these noise distributions were indistinguishable when bias (and/or leakage) was adjusted to give approximately equal mean rate and CV [15]. The examples in Figure 4 show models that simulate the firing normally observed in cat (top row) and firing with a relatively low CV (bottom row). The distributions in Figure 4 are fit with two theoretical curves: the log normal distribution and the diffusion approximation to the random walk with drift. These distributions generally fit better than others that have been used for interval distributions E1.51,including the hyperbolic normal distribution [161, the hyperbolic gamma distribution [25], and the gamma distribution [l, 10,231. The random walk approximation typically fits the simulated data slightly better than the log normal distribution, which is perhaps to be expected because the model (in the absence of leakage) is equivalent to the random walk with drift model. However, the random walk approximation also works well when the integrateand-fire is made leaky. Note that the fits shown for each row of Figure 4 are identical; instead of the best fits to each, the mean best fit is used for all three distributions with a similar CV. Figure 5 shows that the interval distributions generated by the model are appropriate for data. These interval distributions from cat ganglion cells are strikingly similar to the simulations in Figure 4. The log normal distributions and random walk approximations shown on these distributions are identical to those used in Figure 4. In fits to data, the log normal distribution often provides a better fit than the random walk approximation; it also is more forgiving of irregularities such as a possible weak second mode in the interval distribution. Given the ambiguity about the true nature of the distribution, it is probably wiser to use the log normal distribution, for it does not imply particular parameters of the spike-generating process [ 151. While it appears that the interval distributions do not reflect the “signature” of the noise distribution used to generate them, it remains possible that some trace of that distribution can be found in higher moments of the interval distributions. The third moment about the
MQDELING THE VARIABILITY OF FIRfNG RATE
233
234
MICHAEL
Interval
W. LEVINE
(ms)
FIG. 5. Interval distributions from two cat ganglion cells. The smooth curves are identical to those in the corresponding rows of Figure 4. Top: Off-X cell, CV = 0.76 at 21.6 i/s. Bottom: On-X cell, CV = 0.37 at 61.8 i/s. Both recorded from optic tract
161.
mean
is the skew, S:
s=/=(T-?)‘p(+k
(2)
--o
Since P(T) = 0 for T =G0, we may change with some rearranging and substitution, s=
the lower integration
llP3PWd~ 03
1 --
--
(y”
limit, and
3 CV’
Substitution of each of the putative distributions as ~(7) in Eq. (3) gives relationships that can be solved for skew in terms of CV [26]. Specifically, for the random walk approximation (as given in [3]), S random walk= For the log normal
distribution S lognormal
3 cv.
(4)
(as given in [15]), =3cv+cv”.
(5)
MODELING
Finally,
THE VARIABILITY
for the hyperbolic
OF FIRING
gamma
shyperbolic
distribution
RATE
235
(suggested
4cv I- =
in [25]),
(6)
~
l-cv*’
These relationships could allow an additional test of the suggested distributions. Notice that the hyperbolic gamma distribution, which in any case provides a poorer fit than the other two, incorporates the absurdity of S + cc as CV + 1, and S < 0 for CV > 1. Skew is plotted as a function of CV in Figure 6. Various models, using different noise distributions, some with leakage and some without, are shown; for each model, CV (and rate) was adjusted by changing bias. These models produce skews that cluster near and below the 3 CV line predicted for the random walk approximation. There is no obvious discrimination according to leakage or threshold. Data from ganglion cells also fall in a similar band on this plot. The three cells from which data were obtained at the widest range of firing rates (and hence of CV) are shown by solid symbols in Figure 6. Data from other ganglion cells [18] fall in a similar band (not shown on the
0.2
0.4
0.6
0.0
1.0
1.2
FIG. 6. Relationship between skew and CV derived from various integrate-and-fire simulations and data from cat. Simulations: asterisks = Gaussian noise distribution, no leakage; stars = first-order gamma noise, no leakage. Other symbols for uniform noise distributions (square = no leakage; circle = 5 ms time constant, triangles = 2 ms time constant, with two different threshold levels represented by erect or inverted triangles). Solid symbols are data from three cat ganglion cells [6]; solid circles = OffX, solid squares = Off-X, solid diamonds = On-Y.
236
MICHAEL
W. LEVINE
plot). For some cells, the data were noticeably lower than the 3 CV line; these cells display a slight tendency toward bimodal interval distributions. A bimodal distribution can have a relatively larger standard deviation and smaller skew. The models all produce skews comparable to those of the data; however, there is a tendency for the Gaussian noise distribution to produce somewhat higher skews, particularly at high CV. (The Gaussian simulations can produce even higher skews than the ones shown in Figure 6.) Higher skew could be due to the extended lower tail (negative values) of the Gaussian distribution, which both the gamma and uniform distributions lack. This would appear to disqualify Gaussian noise, except that the differences are not significant for the number of simulations and data so far analyzed. If the difference is real, it can easily be ameliorated by postulating a weakly rectifying nonlinearity that attenuates the lower tail of the Gaussian noise. While skew seems not to discriminate among models, it may seem to favor the random walk approximation for fitting interval distributions. However, the relatively low skew (at high CV) is an artifact of the limited data available from a neuron. The theoretical curves were derived by integration to infinity; real data do not include those exquisitely rare very long intervals that are nonetheless inherent in the theoretical curve. Indeed, if an interval in excess of several seconds did occur, the experimenter would assume that some mishap had befallen and discard the run. Because of the cubing operation [see Eq. (211, these rare long intervals have a large effect on skew. Numerical integration of the moments of the log normal distribution shows that for CV near 1 (the fit shown in Figure 2b in [15]), skew is less than 75% of its expected value (for the measured CV) when the area is more than 0.999 (first 325 ms). Thus, the predicted skew for a real interval distribution that is a perfect log normal would be the same as that for the theoretical random walk approximation. This effect is most pronounced for high CV, where high skew is expected, and renders the comparisons in Figure 6 useless. Obviously, examination of higher moments would be futile, for they will be even more disrupted. VARLABILITY
OF RESPONSES
Another way of looking at variability was introduced by workers examining responses of cat cortical cells [5, 28, 291. Cortical cells have little or no maintained discharges, but the responses to repeated stimuli are quite variable. Despite considerable scatter in plots of log variance of rate versus log rate, Tolhurst and his coworkers [5, 28, 291 were able to show that the variance of rate of responses was directly proportional to the mean response amplitude (rate). This was later shown to also be
MODELING
THE VARIABILITY
OF FIRING
RATE
237
true for lateral geniculate neurons (the thalarnic cells that precede cortical cells in the signal pathway from retina to brain) [4]. This result is inconsistent with the nearly inverse square root relationship between CV and mean rate described for cat ganglion cells 114, 231, embodied in Eq. (1). If the firing is a renewal process, as it very nearly is in cat ganglion cells [6], then
where T, is the duration of the samples in which rate is measured [12]. Substitution of Eq. (1) into Eq. (7) shows that variance of rate should be independent of mean rate for any given sample duration [17, 191. It was therefore suggested that the dependence of variance of rate upon mean rate in cortical cells represents a difference between ganglion cells and cortical cells. To test whether there really is a difference between cortical cells and ganglion cells, Levine et al. [191 measured the variance of rate of goldfish ganglion cells in responses to repeated stimuli. They found that variance of rate clearly increased as a function of mean rate, contrary to the expectation from maintained discharges (of cat ganglion cells). Moreover, since ganglion cells support a moderate mean maintained discharge rate, variance of firing was also measured during the phases of the responses that represent decreases in firing rate; variance of rate was a function of firing rate during the response, and not of the strength of the stimulus that elicited that response. There were also some differences from the findings from cortical cells: the variance of rate increased less than in direct proportion to the rate, and the exact nonlinear relationship between variance of rate and rate depended on the length of the periods over which the responses were sampled. In particular, a plot of the logarithm of variance of rate versus the logarithm of mean rate showed that the data were reasonably well fit by a straight line with a fractional slope. The longer the period after each stimulus onset or offset in which responses were measured, the steeper the slope of the best-fit line (see Table 1). This was particularly puzzling: Why would responses measured for longer periods (which thereby include more of the steady maintained-like “plateau”) show more nonlinear interaction with noise than those measured in brief periods (the “purer” responses)? If the variance of rate of maintained discharges is independent of rate [Eq. (7)], the longer durations should show less dependence on mean rate, not more. Since the measurements of variance of rate of ganglion cells were
238
MICHAEL W. LEVINE
made using a different experimental animal, the results did not directly address the question of whether the apparent contradiction between the outcome of the cortical experiments on repeated stimuli and the ganglion cell experiments on maintained discharges was a difference between retinal versus cortical cells or a difference between responses to abrupt stimulation versus maintained discharges. To settle this issue, similar experiments on variability of responses were performed on cat ganglion cells; the results were nearly identical to those from the fish [18] (see Table 1). The only deviation noted was that about one-third of the cells tested showed a distinctly curved relationship on the log variance versus log rate plots. The species did not account for the contradiction; the difference appears to be between responses to sudden transitions in lighting versus maintained discharges in temporally steady conditions. The apparent contradiction between the experiments on variance of rate and those on maintained discharges raised the specter of two kinds of variability that are handled differently in the retina [19]. The noise that affects maintained discharges is added to the rate-setting signal (bias); the noise responsible for the variability of responses to abrupt transitions seems to have a nonlinear relationship with the rate-setting mechanism. How might the retinal network add noise to the steady signal setting maintained rate while multiplying the noise by the transient signals effecting responses ? And, as asked above, why is the multiplication more evident for longer sample durations, which should be more like the maintained discharges? Recent work has shown that these results can be reconciled [17]. As the nonlinearity of the double logarithmic plot derived from steady illumination (Figure 1) illustrates, CV is not actually proportional to the inverse square root of mean rate. Moreover, the firing rate during a response to an abrupt change in illumination of the retina extends through a range of values that may be quite different from the mean rate during the period in which the response is measured. Typically, there is a peak (increase or decrease) immediately upon stimulation, followed by a gradual relaxation to a plateau nearer the maintained firing rate. The longer the duration of the measurement period, the more different the extreme peak will be from the overall mean. The stronger the stimulus, the more extreme the difference will be and the more sharply the peak will decline to the plateau. Levine and Zimmerman [17] proposed a model that, by taking these facts into account, reconciles the two kinds of data. The premise of our model is that the same additive noise affects responses to sudden changes as produces the variability of maintained discharges. The variability is represented by CV, which was assumed to
MODELING
THE VARIABILITY
OF FIRING
RATE
change as a function of rate (F) in the way shown in Figure specifically, a parabolic relationship was assumed:
239
1;
The coefficient of variation is computed from the responses in the response histogram, which is measured in short bins (synchronized to the stimulus and averaged over all the repetitions of the stimulus). This histogram captures the time course of the response, tracing its evolution from peak to plateau. One could not measure CV in such short periods (during a single cycle of which only one or two impulses might occur), but from the mean rate and Eq. (8) one can infer the CV that would be expected. From CV, one can calculate the variance of count expected in that bin; assuming independence of the bins that comprise a response, the variances of count simply sum, and variance of rate can be found by dividing by duration squared. Thus, the variance of rate derived from the mean responses in short time periods can be compared to the observed variance of rate. The model was fit to the data from ganglion cells [18, 191, with the coefficients of the parabolic fit of log CV versus log rate as the only free parameters.3 The model was able to fit the data at least as well as linear regression (despite having fewer free parameters when fitting data from goldfish), even reproducing some of the “scatter” evident in the plots of log variance of rate versus log rate. In those cases in which cat ganglion cells produced curved relationships, the model reproduced the curvature [17]. The mean parabolic relationship from the cat models is plotted as a solid line in Figure 1. It is of the right general form of fit the observed maintained discharges, which are well fit by the leaky integrate-and-fire model. The vertical error is not statistically significant and is necessarily linked to the actual CV of the sample (see footnote 3). Remember that the maintained discharge data and the data upon which the response
3Equation (8) has three coefficients, but only two are used as parameters of the model. A is found from X, and .r2 given the mean rate and CV of the observed maintained discharge of the cell. The third model parameter is a multiplier used in converting CV to variance of rate; this fractional factor (very nearly unity) accounts for the slight temporal structure of the firing, which is not exactly a renewal process 16, 12, 131. The value of the third parameter is consistent with the values in rows 8 and 9 of Table 1. Since multiplication of CV is equivalent to adding a constant to log CV, this is equivalent to allowing adjustment of A. However, the values of A used for the smooth curve in Figure 1 remain set to reproduce the observed CV of each cell.
240
MICHAEL
W. LEVINE
variability model were based were drawn from experiments, in laboratories hemispheres apart, thetics.
two different sets of using different anes-
CONCLUSION Our basic model suggests that the variable discharges of a ganglion cell are generated by adding noise to the deterministic signal representing visual information, and using this combined signal as the input to a leaky integrate-and-fire mechanism. This model accounts well for the variability of maintained discharges, both in terms of the distribution between successive neural impulses and in the way in which variability changes as mean rate shifts in response to steady stimulation. It also can reproduce the higher order statistics and temporal dependencies in the firing [6, 13, 141. It is also consistent with the variability of firing in response to transient stimulation. data, it must be Although the model accounts for “well-behaved” noted that interval distributions from ganglion cells are not always smoothly unimodal. A minority of ganglion cells produce multimodal interval distributions, and these cannot always be dismissed as pathological (as some workers do). The models considered here make no attempt to explain multimodal distributions, although a modification in which the threshold changes after each impulse can produce such an effect [14]. Alternatively, multimodal distributions could be due to a “multiplexing” among different noise or signal sources [16] or to a second mechanism by which impulses may be produced according to the derivative of the membrane potential [ll, 20, 211. These complications have not been considered here. Finally, it is important to remember that even if an integrate-and-fire model is satisfactory, it is not unique. The trade-offs among parameters allow the same data to be fit with a large number of models; intracellular measurements will be required to determine the actual parameters of ganglion cells. I am gratefil to the many collaborators over the years who have offered advice and suggestions as this work progressed; major contributions were made by Violeta Cam’&-Carire, Jeremy Shefner, and especially Roger Zimmerman. I am particularly grateful to those collaborators who have allowed me to share their facilities and harness their expertise to collect the data reported here, speci@ally Laura Frishman and Christina Enroth-Cugell at Northwestern University and Brian Cleland at the University of Sydney. I also wish to express my appreciation to Grace Yang and Charles Smith, organizers of this conference, for inviting me to participate.
MODELING
THE VARIABILITY
OF FIRING
RATE
241
REFERENCES
4
10 11 12 13
H. B. Barlow and W. R. Levick, Changes in the maintained discharge with adaptation level in the cat retina J. Physiol. 202:699-718 (1969). H. B. Barlow, W. R. Levick, and M. Yoon, Responses to single quanta of light in retinal ganglion cells of the cat, I/is&n Res. ll(Supp1. 3):87-101 (1971). D. Berger, K. Pribram, H. Wild, and C. Bridges, An analysis of neural spike-train distributions: determinants of the response of visual cortex neurons to changes in orientation and spatial frequency, Exp. Brain. Res. 80:129-134 (1990). V. Carrion-Carire, B. G. Cleland, A. W. Freeman, M. W. Levine, and R. P. Zimmerman, Variability of responses in retinal and geniculate cells, Int,est. Ophthalmol. Vis. Sci. 29(Suppl):295 (1988) (Abstract). A. F. Dean, The variability of discharge of simple cells in the cat striate cortex, Exp. Bruin Res. 44:437-440 (1981). L. J. Frishman and M. W. Levine, Statistics of the maintained discharge of cat retinal ganglion cells, J. Physiol. 339:475-494 (1983). G. L. Gerstein and B. Mandelbrot, Random walk models for the spike activity of a single neuron, Biophys. J. 4:41-68 (1964). P. I. M. Johannesma, Diffusion models for the stochastic activity of neurons, in Neural Networks, E. R. Caianiello, Ed., Springer-Verlag, Berlin, 1968, pp. 116-144. B. W. Knight, Dynamics of encoding in a population of neurons, J. Gen. Physiol. 591734-766 (1972). S. W. Kuffler, R. FitzHugh,
and H. B. Barlow, Maintained activity in the cat’s retina in light and darkness, J. Gen. Physiol. 40:683-702 (1957). M. J. M. Lankheet, J. Molenaar, and W. A. van de Grind, The spike generating mechanism of cat retinal ganglion cells, Vision Res. 29:.505-517 (1989). M. W. Levine, Firing rates of a retinal neuron are not predictable from interspike interval statistics, Biophys. /. 30:9-26 (1980). M. W. Levine, Retinal processing of intrinsic and extrinsic noise, J. Neurophysiol. 48:992-1010
(1982).
M. W. Levine, Variability in the maintained discharges of retinal ganglion cells, J. Opt. Sot. Am. A 4:2308-2320 (1987). I5 M. W. Levine, The distribution of the intervals between neural impulses in the maintained discharges of retinal ganglion cells, Biol. Cybem. 65:459-467 (1991). 16 M. W. Levine and J. M. Shefner, A model for the variability of interspike intervals during sustained firing of a retinal neuron, Biophys. J. 19:241-252 14
(1977).
17 18 19 20
M. W. Levine and R. P. Zimmerman, A model for the variability of maintained discharges and responses to flashes of light, Biol. Cybem. 65~469-477 (1991). M. W. Levine, B. G. Cleland, and R. P. Zimmerman, Variability of responses of cat retinal ganglion cells, Visual Neurosci. 8:277-279 (1992). M. W. Levine, R. P. Zimmerman, and V. Carrion-Carire, Variability in responses of retinal ganglion cells, /. Opt. Sot. Am. A 5:593-597 (1988). R. Llinas and H. Jahnsen, Electrophysiology of mammalian thalamic neurones in vitro, Nature 297:406-408
21
(1982).
F.-S. Lo and S. M. Sherman, In vivo recording of postsynaptic potentials and low threshold spikes in W cells of the cat’s lateral geniculate nucleus, Exp. Bruin Res. 81:438-442 (1990).
242 22 23 24 25 26 27 28
29
MICHAEL
W. LEVINE
J. Pernier and P. Gerin, Temporal pattern analysis of spontaneous unit activity in the neocortex, Biol. Cybem. 18:123-136 (19751. J. G. Robson and J. B. Troy, The nature of the noise in cat retinal ganglion cells, J. Opt. Sot. Am. A 4:2301-2307 (1987). R. W. Rodieck, Maintained activity of cat retinal ganglion cells, J. Neurophysiol. 30:1043-1071 (1967). C. E. Smith, A comment on a retinal neuron model, Biophys. /. 25:385-386 (1979). C. E. Smith, personal communication (1991). R. B. Stein, A theoretical analysis of neuronal variability, Biophys. J. 5:173-194 (1965). D. J. Tolhurst, J. A. Movshon, and A. F. Dean, The statistical reliability of signals in single neurons in cat and monkey visual cortex, Ksion Res. 23:775-785 (1983). D. J. Tolhurst, J. A. Movshon, and I. D. Thompson, The dependence of response amplitude and variance of cat visual cortical neurones on stimulus contrast, Exp. Bruin Res. 41:414-419 (1981).
