Previous: Modeling Neurophysiological Data Up: Modeling Neurophysiological Data Next: Saturation levels revisited

Reconstruction of 3D response functions

The data in [De Valois et al. 1966] consist of the averaged responses of six types of cells in the Macaque LGN to monochromatic stimulus light, sampled at three radiance levels and 12 to 13 wavelengths for each cell type. The cell types are designated according to the approximate color of the light that results in a maximum excitatory (+) or inhibitory () response: the four opponent types , , and , and the two non-opponent types referred to as excitators and inhibitors, which I will refer to as and , respectively (the symbols are for red, for green, for yellow, for blue, and for any light). Figure shows the six data sets.

The inhibitory (negative) phases of the responses are generally much smaller than the excitatory (positive) phases, because they are relative to the spontaneous firing rate, which is typically 5-10 spikes/s, and negative firing rates are obviously not possible. The columns of Figure show pairs of cell types that are complementary in nature, with responses that are a kind of ``mirror image'' of each other (when disregarding the absolute sizes of the phases): in a wavelength range where one type is excited, the other is inhibited (and the other way around), and the cross-over wavelengths between excitation and inhibition (if any) are virtually identical. The use of these mirror-imaged channels is common in the vertebrate visual system, and may be a way to encode a larger dynamic range than would be possible with a single channel [Slaughter 1990].

To reconstruct the 3D response functions from this relatively sparse data set, I started by fitting a model to the data in the radiance domain, followed by interpolation in the wavelength domain. The reason for starting in the radiance domain is that a lot is known about the response characteristics of photoreceptors and other cells in the visual system with respect to changing stimulus intensity, and this knowledge provides important constraints on the model. Based on the literature (e.g., [Frumkes 1990][Slaughter 1990][Leibovic 1990a]), I have assumed a sigmoid-shaped function as the basic response model in the radiance domain:

where represents radiance, the half-response radiance, and is a parameter that determines the first derivative around the half-response point, sometimes referred to as the temperature parameter. Figure shows a plot of this function for and .

In practice I use a scaled version of the sigmoid function:

where c is a scaling factor that is usually 1, but takes on different values on some occasions (see below).

This type of function is particularly appropriate for modeling, e.g., the operating curve of photoreceptors that shifts along an absolute radiance axis as a result of adaptation, as background radiance increases or decreases [Leibovic 1990a]. The visual system's response is always relative to its adaptation state, and if we assume a fixed adaptation state, a function like the one in Figure can be used to model the response. The parameter can be interpreted as a function of the adaptation state. Near , the sigmoid function is almost linear, but moving away from towards and it approaches the limits and , respectively. This can be interpreted as approaching threshold response and saturation response, respectively. The sigmoid function is also commonly used as a node activation function in artificial neural networks research (e.g., [Rumelhart et al. 1986]).

To model the response of an LGN cell (type) as a function of radiance (for a fixed wavelength) with a sigmoid function, we need to interpret the data points in terms of and , and then determine appropriate values of and that minimize the error of fit. To make the computations easier, I will use the interval as the ``useful'' range of corresponding to a value going from near-zero (threshold) to near-one (saturation). The radiances in [De Valois et al. 1966] are specified as log units below some maximum used radiance. Log transformations of physical intensity are generally used in psychophysics to express perceived intensities, since they roughly correspond to a perceptually linear scale, over a significant part of the useful range. Since these radiances are relative, we need to decide where they should be situated along the radiance axis. I made the following assumptions:

1. The dynamic range of the responses from threshold to saturation is 3.5 log units of input radiance. This accords well with estimates of 3 log units for the ``useful'' (near-linear) range given in [Leibovic 1990b] for photoreceptors, extended slightly to allow for response threshold and saturation. This value works well for fitting the sigmoid intensity response model to the data, as explained below. Since we know the difference between the highest and lowest radiance level used in the experiments, and all measurements are made at a constant adaptation state, we can derive a rough estimate of the maximum level used relative to the background light level, as

   

   where is the maximum level used in log units relative to the background light level, is the highest and is the lowest stimulus radiance used relative to the maximum available, is the average firing rate across the spectrum at the highest and at the lowest radiance used. I used the average responses to minimize measurement errors. Of course this equation assumes linear response as a function of log radiance, which we know not to be the case, but to be approximately true around the half-response radiance. Table summarizes the results.

   

   In some cases (viz. +R-G, +B-Y, and especially -L), the estimated maximum is smaller than the maximum attenuation used, which would imply light levels below that of the background adaptation light. The article is not clear on whether or not this was the case, and I have not been able to determine this any other way. The -L data are all inhibitory, i.e., negative with respect to the spontaneous firing rate, and thus has a small range. There may hence be a considerable measurement error involved for this data set. In any case, it is safe to assume that the total useful range of input radiances must be at least as large as the average of the values as given above. In fact it must be at least as large as the maximum value of (the maximum attenuation used), viz. 1.7 log units, since the cells are still responsive at that level. The latter argument would give us a range of about 2 log units. Since the measurements were most likely done at radiance levels well below saturation, the total range must be larger than this, which is consistent with the figure of 3.5 log units mentioned above.

2. The preceding argument gives us a radiance range of about 2 log units for the experimental data, with a maximum of 3.5. The radiance of the brightest stimuli is therefore about 1.5 log units below saturation level. The radiance levels for each of the data sets are shown in Figure (p. ), normalized to .

3. The responses of all six cell types will eventually saturate at approximately the same level if the stimulus radiance is turned up high enough, since this is more a function of the biochemistry of the cell than of the absolute stimulus intensity. To scale the responses of the cell types, I used a global maximum (saturation level) of 150 spikes/sec above the spontaneous firing rate. This is somewhat lower than typical LGN cell saturation levels of 2-300 spikes/sec, but the spontaneous rate in these data of 5-10 spikes/sec is also lower than the typical rates of about 30. The fact that the experimental subjects were anesthetized may have something to do with this general lowering of the firing rates. The maximum occuring firing rate in the data of about 50 spikes/sec for the +Y-B cell type at about 610 nm for an estimated stimulus radiance of about 2 log units would support a saturation level of about 150 spikes/sec for a total of 3.5 log units, taking the non-linearity of the response at high radiances into account.

4. Since I measure responses relative to the spontaneous firing rate, I also need to determine a maximum negative response. These responses represent inhibition, and the absolute firing rate can obviously not drop below 0. Since the spontaneous rates in De Valois' data does not exceed 11 spikes/sec, I've set the negative maximum to -12.

Although the estimates of the total useful stimulus intensity range and saturation levels are somewhat speculative (since there is not enough data to support better estimates), they are generally not in disagreement with known characteristics of LGN cell responses. In addition, the exact numbers chosen are not all that important, since using different values amounts basically to a linear scaling of the (non-linear) response functions, as I will show below.

With these assumptions in place, we can now fit the sigmoid models to the data. I manually extrapolated the data in the lower and upper wavelength ranges, so each constant-radiance data set starts and ends with a zero point. In addition, I added a zero point to each constant-wavelength set of data points, to constrain the sigmoid fit better. I used the following algorithm to do the fit for each constant-wavelength data set:

 
1.  let data = sign(data)*data 
1.1 c = {Abs[NegMax/PosMax], sign < 0 
        { 1.0 , otherwise 
2.  compute least-squares linear fit of data as function of r: a=mr+b 
3.  let h1 = solution(mr+b=0.5 c) solving for r 
4.  let t1 = solution(S'(h1,h1,t,c)=m) solving for t 
5. minimize(SE2(r,h1,t1,c)) solving for h1 and t1 
6.  return sign(data)*S(r,h1,t1,c)

Since the sigmoid is really a function of three variables (with radiance , half-response radiance , and temperature , cf. equation p. ), we need to optimize the fit with respect to one, , by varying the other two, and . The data form a list of ordered pairs for a fixed wavelength and cell type, including the added point. Step 1 makes sure that all response values are positive. In principle, a data set may contain both negative and positive response values, especially around the cross-over point between inhibition and excitation, although in practice this is rare. This step is more useful for the inhibitory phases of the response, which are represented as negative values since I measure responses relative to the spontaneous rate. I define the sign of the data set as the sign of the sum of the response values in the set:

where is the data set, is the function that selects the second element from an ordered pair, is the relative radiance, and is the response (activation). The sigmoid function has to be either entirely positive or entirely negative, and this step chooses between those. Note that this also turns a set of data points whose 's sum to zero into all-zero 's. This can only happen around the cross-over point between inhibition and excitation, and is not wrong, even desirable, because varying the radiance alone should not affect the sign of the response, so the data are probably unreliable in that case. Step 1.1 chooses the scaling factor for the sigmoid function, 1 for positive data and the ratio of the maximum negative to the maximum positive response for negative data. To find initial values for and , we first compute a linear fit to the data set in step 2, using Mathematica's least-squares Fit function. Solving for a linear equation value of gives us an initial estimate for in step 3 (recall that is the half-response radiance, and the sigmoid function values are normalized to ). Step 4 gives us an estimate for such that the first derivative of the sigmoid function at is equal to the slope of the linear equation. This is a reasonable estimate, since the sigmoid function is near-linear around the half response point, and the data points cover approximately the first half of the total radiance range. Step 5 does the actual fit, using Mathematica's FindMinimum function (which finds local minima using a steepest gradient descent algorithm) on SE2, which is the summed squared error of fit over the data set. The last step returns the sigmoid function model for the data set, restoring the sign. The algorithm works well in practice. Figure illustrates the procedure, and Figure shows some examples of the computed sigmoid fits.

Actually, steps 2 through 4 can be omitted and the minimization started with initial estimates of for both and , but in that case it may take longer to converge. A constraint is imposed on the fit; if it is violated, a different iterative error minimization procedure is used that keeps . This constraint serves to prevent slopes of the sigmoid function that are too steep to serve as realistic models of neural activation, but in practice this is seldom necessary.

After computing an array of sigmoid fit functions (indexed by wavelength) for each sampled wavelength and each of the six cell types, we can compute the response of each type to an arbitrary wavelength and radiance by interpolating between the values given by the sigmoid functions of radiance that are closest to the given wavelength:

where is the index into the array of sigmoid functions that we are looking for, is the wavelength for which we want to compute a response, and is the -th wavelength for which we have computed a sigmoid function. Note that before this interpolation is applied, there is no smoothness constraint imposed on the parameters of the sigmoid fit from one sampled wavelength to the next, so in principle there could be serious discontinuities. In practice one would expect the parameters to change smoothly, and it turns out they do, to a considerable extent at least. I used Mathematica's standard cubic spline Interpolation function on the 7 sigmoid function values returned by the series of functions centered around , adding additional zero points above and below the sampled wavelength range (and specifying 0 first and second derivatives for the very lowest and highest, to constrain the interpolation better):

where is the computed response (activation), computes a cubic spline interpolation function for the points in list and evaluates it at , is the -th radiance level and the corresponding sigmoid function we have previously computed, and is the relative radiance and the wavelength for which we want to compute the response. The resulting functions are shown in Figure for the same relative radiances as in Figure (p. ), and additionally for maximum relative radiance .

We can see that the functions follow the data sets closely at the sampled wavelengths and radiances, but they are continuous and extend beyond the sampled data in both domains. We can see some irregularities in the maximum radiance level responses, which results from extrapolating the data considerably, combined with the absence of smoothness constraints as explained above. These irregularities will largely disappear in the next processing step, as explained in Section (p. ). It is interesting to note the clear saturation effect in most inhibitory response phases, and in the excitatory phases of the and functions, due to the use of the sigmoid model. If we continued to compute function values at ever higher radiances, saturation would eventually also set in in the other functions. It is not clear whether the radiance range should be extended further to allow all functions to reach saturation levels. Alternatively, the radiance axis could be normalized for each function individually, or perhaps for opponent pairs of functions, but that is not in agreement with the assumptions I outlined above. Saturation effects are clearly in agreement with neurophysiological data, but to my knowledge there is no other color model available that incorporates them.