Connecting psychophysical performance to neuronal response properties II: Contrast decoding and detection

Abstract

The purpose of this article is to provide mathematical insights into the results of some Monte Carlo simulations published by Tolhurst and colleagues (Clatworthy, Chirimuuta, Lauritzen, & Tolhurst, 2003; Chirimuuta & Tolhurst, 2005a). In these simulations, the contrast of a visual stimulus was encoded by a model spiking neuron or a set of such neurons. The mean spike count of each neuron was given by a sigmoidal function of contrast, the Naka-Rushton function. The actual number of spikes generated on each trial was determined by a doubly stochastic Poisson process. The spike counts were decoded using a Bayesian decoder to give an estimate of the stimulus contrast. Tolhurst and colleagues used the estimated contrast values to assess the model's performance in a number of ways, and they uncovered several relationships between properties of the neurons and characteristics of performance. Although this work made a substantial contribution to our understanding of the links between physiology and perceptual performance, the Monte Carlo simulations provided little insight into why the obtained patterns of results arose or how general they are. We overcame these problems by deriving equations that predict the model's performance. We derived an approximation of the model's decoding precision using Fisher information. We also analyzed the model's contrast detection performance and discovered a previously unknown theoretical connection between the Naka-Rushton contrast-response function and the Weibull psychometric function. Our equations give many insights into the theoretical relationships between physiology and perceptual performance reported by Tolhurst and colleagues, explaining how they arise and how they generalize across the neuronal parameter space.

Introduction

This work was inspired by a series of articles from Tolhurst and colleagues that described Monte Carlo simulations of contrast coding in the visual cortex (Chirimuuta, Clatworthy, & Tolhurst, 2003; Clatworthy, Chirimuuta, Lauritzen, & Tolhurst, 2003; Chirimuuta & Tolhurst, 2005a, 2005b). In these simulations, the contrast of a visual stimulus was encoded by a model spiking neuron, or a set of such neurons, and the spike counts were decoded using a Bayesian decoder to give an estimate of the stimulus contrast. The contrast estimates were then used to assess the model's decoding precision or to generate decisions in simulations of psychophysical tasks. 

Although these studies made a very useful contribution to our understanding of the relationship between physiology and perceptual performance, the Monte Carlo simulations provided little insight into why the obtained patterns of results arose or how general they are. It was not clear whether the findings applied to just the specific sets of model parameters that were used in the simulations or to any parameter values. The goal of the present study was to derive equations that predict the performance of this kind of model from its parameters in order to better understand the relationships between physiology and perceptual performance and to allow the model's performance to be calculated quickly and easily, thus removing the need for Monte Carlo simulations, which are difficult to program and time consuming to run. 

Our article focuses on two key sets of results from the work of Tolhurst and colleagues. The first focus is the relationships between decoding precision and neuronal properties discovered by Clatworthy et al. (2003); we find that these relationships can be explained by approximating the decoding precision using Fisher information. The second focus is the model's psychometric function for detection, studied by Chirimuuta and Tolhurst (2005a); we derive an expression that approximates this function and sheds light onto one of Chirimuuta and Tolhurst's (2005a) results. 

Some of the derivations and technical details are included in appendices in the supplementary information. Each supplementary appendix is labeled with a letter. Equations and figures in an appendix are labeled with the appendix's letter, followed by a dot, followed by the equation or figure number within that appendix. Supplementary Appendix A provides a list of the main symbols used in the text and their meanings. 

The contrast coding model

The contrast coding model used by Tolhurst and colleagues consisted of an encoding stage, in which the stimulus was encoded as neuronal spike counts, and a decoding stage, which used Bayesian inference to estimate the stimulus contrast from the spike counts. The encoding stage was based on standard physiological models of the neuronal response, while the decoding stage was simply a Bayesian decoder, with no proposed physiological implementation. 

The model has three elements: (a) the form of the tuning function, which specifies the mean spike count of each neuron for a given stimulus; (b) the random process that generates spikes at the given rate; and (c) the method of decoding the population response to give an estimate of the stimulus. We consider these three elements in turn. 

We use the same basic terminology as May and Solomon (2015): We represent random variables using uppercase letters and their values on particular trials¹ using corresponding lowercase letters. X and x represent the stimulus level, R and r(x) represent the mean spike count, N and n represent the actual spike count, and N and n are vectors holding the spike counts of all the neurons in the population. These variables are explained in more detail in our companion article (see the first paragraph of the section titled “The sensory coding model” in May & Solomon, 2015). 

The tuning function

The tuning function, r(x), specifies the mean spike count for stimulus x. For visual stimulus contrast, the tuning function is called the contrast-response function. Tolhurst and colleagues modeled this function using the Naka-Rushton function (Naka & Rushton, 1966; Albrecht & Hamilton, 1982). For contrast, c, in linear (e.g., Michelson) units, the Naka-Rushton function has the following form:  

 

If we measure contrast in log units, x = log_bc, then the Naka-Rushton is given by  

 

where  

 

See our companion article (May & Solomon, 2015) for a description of all the parameters and for plots of the Naka-Rushton function. On the log contrast scale, the gradient of the Naka-Rushton function peaks at a log contrast of x = z. The term c_1/2, called the semisaturation contrast, is the contrast for which the mean response exceeds r₀ by r_max/2. Many authors use c₅₀ to represent this contrast, but we use the less common term c_1/2 (except when quoting other authors) because this form is easier to extend to other fractions of r_max, such as c_1/3 (the contrast for which the mean response exceeds r₀ by r_max/3, which we show to be the contrast that is coded most accurately by a model neuron with a Naka-Rushton contrast-response function with r₀ = 0).  

In this article, whenever we use the term contrast without specifying the units, we mean log Michelson contrast. To be compatible with Clatworthy et al. (2003) and Chirimuuta et al. (2003), we always used log to base 10 in our modeling (i.e., x = log₁₀ c and z = log₁₀ c_1/2); however, our equations are derived for the general case of any base, b. 

Note that the units of the Naka-Rushton function's output are often taken to be spikes per second, but we find it more convenient to use units of spikes without making assumptions about the time interval over which the neuron's output is integrated. 

The random process for spike generation

The tuning function gives the mean spike count, and we now turn to the stochastic process that generates spikes at the specified rate. The Poisson process is sometimes used as a stochastic spiking model because of its considerable mathematical convenience (e.g., see Dayan & Abbott, 2001, chapter 1). The Poisson distribution is defined as  

 

In Equation 4, we use a right subscript on P to indicate the type of stochastic process being used.  

For the Poisson process, the variance in the number of spikes for repeated trials is equal to the mean, r(x), always giving a Fano factor (ratio of variance to mean) of 1. In the visual cortex, the Fano factor is usually greater than 1 (Dean, 1981b; Tolhurst, Movshon, & Thompson, 1981; Tolhurst, Movshon, & Dean, 1983; Bradley, Skottun, Ohzawa, Sclar, & Freeman, 1987; Skottun, Bradley, Sclar, Ohzawa, & Freeman, 1987; Scobey & Gabor, 1989; Vogels, Spileers, & Orban, 1989; Snowden, Treue, & Andersen, 1992; Britten, Shadlen, Newsome, & Movshon, 1993; Softky & Koch, 1993; Swindale & Mitchell, 1994; Geisler & Albrecht, 1997; Bair & O'Keefe, 1998; Buracas, Zador, DeWeese, & Albright, 1998; McAdams & Maunsell, 1999; Oram, Weiner, Lestienne, & Richmond, 1999; Durant, Clifford, Crowder, Price, & Ibbotson, 2007). To get a Fano factor greater than 1, Tolhurst and colleagues (Chirimuuta et al., 2003; Clatworthy et al., 2003; Chirimuuta & Tolhurst, 2005a, 2005b) used a doubly stochastic Poisson process, which we refer to as the Tolhurst process. This process is a Poisson process in which the mean is itself a random variable sampled from a Poisson process with mean r(x):  

 

For this process, the mean spike count is r(x) and the variance is 2r(x), giving a Fano factor of 2. The infinite series in Equation 5 is difficult to handle, so in Supplementary Appendix B we derive a finite series expansion of the Tolhurst process that is more useful.  

Since the aim of this article is to explain the modeling results of Tolhurst and colleagues, all the Monte Carlo simulations in this article use the Tolhurst process to generate spikes. However, to make our theoretical results more general, we consider two other Poisson-based spiking processes. 

The first additional spiking process is the generalized Poisson distribution, devised by Consul and Jain (1973). Sakata and Harris (2009) used this process to model neuronal spiking distributions. This process is more flexible than the Tolhurst process because the Fano factor, F, is set as a parameter and can take any value greater than or equal to 1. The Consul-Jain distribution takes the following form:  

 

where  

 

Like the Poisson and Tolhurst processes, the Consul-Jain process generates only nonnegative numbers of spikes, and the spike count variance is proportional to the mean. When F = 1, the Consul-Jain process reduces to the ordinary Poisson process, given in Equation 4. However, the Consul-Jain process with F = 2 is not identical to the Tolhurst process, even though the Fano factor is the same.  

The second additional spiking process that we consider is the doubly stochastic Poisson process recently proposed by Goris, Movshon, and Simoncelli (2014). We refer to this as the Goris process. It is a Poisson process in which the mean is modulated by a multiplicative gain mechanism. The gain value on each stimulus presentation is a gamma-distributed random variable. The Goris process is described in detail in our companion article (May & Solomon, 2015). The gain fluctuations result in a Fano factor that increases with the firing rate, which is more physiologically plausible than the fixed Fano factor produced by the Tolhurst and Consul-Jain processes. For mathematical simplicity, we assume that each neuron in the population has the same gain signal and that each neuron's Poisson spiking process is independent (see May & Solomon, 2015, for justification of these restrictions). The shared gain signal causes the neuronal responses to be correlated, with a realistic correlation structure. However, because each neuron's Poisson spiking process is independent, a decoder that knows the gain signal on each stimulus presentation can express the spiking distributions as independent Poisson distributions. Therefore, the neurons can be decoded as if they were statistically independent (see May & Solomon, 2015, for a more detailed explanation of this). 

Bayesian population decoding

All the results that we report from Tolhurst and colleagues were obtained using maximum-likelihood decoding. The estimated stimulus level is the one that had the highest probability of giving rise to the obtained set of spike counts; that is, it is the value of x that maximizes the likelihood, P(N = n|X = x). In Tolhurst and colleagues' model, the neurons were statistically independent; in addition, if we use May and Solomon's (2015) parameterization of the Goris process, the neurons are implicitly decorrelated if the decoder knows the gain signal. For statistically independent neurons, the population likelihood is then given by the product of the likelihoods of the individual neurons:  

 

 

 

K is the number of neurons. N_j is a random variable representing the spike count of neuron j, and n_j is its value. R_j is a random variable representing the mean spike count for neuron j, and r_j(x) is its value. The second equality (Equation 9) arises because each r_j(x) is a deterministic function of the stimulus value, x. To evaluate the probabilities in Equation 9, we use the appropriate expression, depending on which spiking process we are using. For the Consul-Jain process, we use Equation 6. For the Goris process, assuming the decoder knows the gain value, we use the gain-modulated Poisson distribution (May & Solomon, 2015; Equation 10). All the modeling in this article used the Tolhurst process, which is defined in Equation 5.  

Measuring decoding precision

Definition of decoding precision

Decoding precision is usually taken to be the reciprocal of the variance of the estimated stimulus value; this is how it is defined in our companion article (May & Solomon, 2015). However, in their Monte Carlo simulations, Tolhurst and colleagues (Chirimuuta et al., 2003; Clatworthy et al., 2003) used a slightly different measure of decoding performance. They called it accuracy, defined as the reciprocal of the mean squared difference between the individual estimate and the true stimulus value:  

 

where T is the number of trials (= 10,000), x̂ is the estimated log contrast estimate (note that we use an uppercase x̂ to represent the random variable and a lowercase x̂ to represent its value on a particular trial), and the denominator is the sum over all trials. For the models that we consider in this study, the log contrast estimate is largely unbiased (except at very low performance levels), so we have mean[x̂] ≈ x. In this case, the accuracy score in Equation 10 is essentially the same as the precision. For consistency with Tolhurst and colleagues, we used their measure of decoding accuracy (Equation 10) when analyzing our Monte Carlo simulations. However, we refer to it as precision (except when directly quoting Tolhurst and colleagues) because our analytical approximations of it are formally measures of precision and because the term accuracy is often used to mean the inverse of bias (e.g., Smith, 1999, chapter 2). We found that, except in degenerate conditions where the model's performance level was very low, it made a negligible difference whether we plotted Tolhursts's accuracy score or true precision. Note that in the Monte Carlo simulations in our companion article (May & Solomon, 2015), we calculated true decoding precision (i.e., reciprocal of the variance of the estimated stimulus value), not the accuracy score defined in Equation 10.  

Approximating decoding precision using Fisher information

For an unbiased maximum-likelihood decoder, the precision cannot exceed a quantity called the Fisher information. This limit is known as the Cramér-Rao bound (Rao, 1945; Cramér, 1946; see Dayan & Abbott, 2001, pp. 120–121). For populations that give sufficiently large numbers of informative spikes, the precision of a maximum-likelihood decoder is very close to the Fisher information; when the tuning function is a sigmoid function like the Naka-Rushton function, this applies even if the population consists of a single neuron. This means that we can use the Fisher information as an analytical approximation of the decoding precision, thereby giving us a real insight into why Tolhurst and colleagues' results occurred and how general they are. 

The Fisher information depends on the tuning function and the spiking process. Unfortunately, for the Tolhurst and Consul-Jain spiking processes, the derivation of an exact formula for the Fisher information turned out to be intractable. In Supplementary Appendix E, we derive approximations of the Fisher information for each of these spiking processes: τ_Tolhurst(x) and τ_C-J(x). As long as the mean spike count of the most informative neurons is not too low, both τ_Tolhurst(x) and τ_C-J(x) are very close to the following general approximation of the decoding precision, τ̃(x):  

 

where Display Formula (x) is the first derivative of neuron j's tuning function with respect to x. To parameterize τ̃(x) so that it approximates τ_Tolhurst(x) or τ_C-J(x), v should be equal to the Fano factor; the Fano factor is fixed at 2 for the Tolhurst process and can take any value greater than or equal to 1 for the Consul-Jain process. The tilde (˜) above τ indicates that this estimate of decoding precision is based on an approximation of the Fisher information that is not always accurate. With a Fano factor of 1, the Consul-Jain process is the ordinary Poisson, in which case τ̃(x) with v = 1 is exactly equal to the Fisher information (see Dayan & Abbott, 2001, chapter 3). For the Tolhurst process, and the Consul-Jain process with F ≠ 1, τ̃(x) with v equal to the Fano factor is an approximation of the Fisher information. For the Goris process, May and Solomon (2015) showed that an appropriate estimate of the decoding precision is given exactly by Equation 11 with v = 1/(1 − Display Formula ), where Display Formula is the variance of the gain signal (see equation 26 and relation 31 of May & Solomon, 2015). In this case, v is not the Fano factor—the Fano factor for Goris et al.'s (2014) spiking process is variable and depends on the mean spike count. When the gain is known by the decoder, the Fisher information of the Goris process varies from trial to trial due to the fluctuating gain, and τ̃(x) with v = 1/(1 − Display Formula ) gives the modal value of the Fisher information exactly.  

In summary, Equation 11 gives a general approximation of the decoding precision for a variety of Poisson-based spiking distributions. For two processes (Poisson and Goris) it is based on an exact expression for the Fisher information; for the other two processes (Tolhurst and Consul-Jain) it is based on an approximation of the Fisher information that is accurate as long as the mean spike count of the most informative neurons is not too low. 

If we expand r_j(x) and Display Formula (x) in Equation 11 using the Naka-Rushton function (Equation 2), then we have  

 

Tolhurst and colleagues always used r₀ = 0 in their modeling. Since we are focusing on their modeling results, we consider only the case of r₀ = 0 in this article. In Equation 12, each neuronal parameter can vary from neuron to neuron, so, strictly speaking, each parameter should be indexed by the neuron number, j, but we omit these indices to reduce notational clutter.  

In the next section, many of the populations of neurons that we analyze consist of either just a single neuron or a set of identical neurons. In this case, Equation 12 reduces to  

 

with K = 1 in the case of a single neuron.  

Explaining Clatworthy et al.'s contrast decoding results

Clatworthy et al. (2003) applied their Bayesian contrast decoding algorithm to single Tolhurst-spiking neurons, or populations of them, and discovered several relationships between the neuronal parameters and the peak of decoding precision along the contrast axis. 

1.  
   “Increasing R_max increases the contrast identification accuracy of single neurons at all contrasts, most obviously the peak accuracy, without changing the contrast at which accuracy is a maximum” (Clatworthy et al., 2003, p. 1991; note that they use R_max where we use r_max). This can be summarized by saying that there is an approximately multiplicative effect of r_max on decoding precision.
2.  
   “The position of the accuracy peak along the contrast axis is consistently close to but, interestingly, slightly below the neuron's c₅₀” (Clatworthy et al., 2003, p. 1989).
3.  
   “The peak value of accuracy is independent of c₅₀” (Clatworthy et al., 2003, p. 1989).
4.  
   “The relationship between the maximum accuracy and q is a steep straight line on log-log coordinates” (Clatworthy et al., 2003, p. 1989).
5.  
   “To change the maximum accuracy . . . requires only a change in the product of R_max and number of neurons, i.e., the total number of action potentials generated on average . . . ; for a given accuracy, there is a simple trade-off between the number of neurons and the response amplitude of individual neurons” (Clatworthy et al., 2003, p. 1990).

The numerical nature of Clatworthy et al.'s method means that it gave little insight into why these relationships arose or how generally they apply. In this section we explain Clatworthy et al.'s findings by deriving equations that explain corresponding findings for the analytical approximation of decoding precision, τ̃(x), given in Equation 13. 

First, we examined whether τ̃(x) was sufficiently close to the true decoding precision to make this approach valid. This is important because the Fisher information, on which τ̃(x) is based, is only an upper bound on the decoding precision and can far exceed the true decoding precision for small population sizes (Xie, 2002). In one of their investigations (shown in their figure 5A), Clatworthy et al. (2003) examined the effect of r_max on decoding precision for a single Tolhurst-spiking neuron with Naka-Rushton exponent (q) equal to 2, c_1/2 = 0.1, and r₀ = 0. r_max took values of 5, 20, 50, 100, or 180 spikes. We replicated their methods (see Supplementary Appendix G for details) and obtained contrast decoding precision scores (calculated using Equation 10) that were essentially identical to those that we read off from their figure 5A. The small differences were almost certainly attributable to the stochastic nature of the simulations and small inaccuracies in our transcription of the data from Clatworthy et al.'s figure. The symbols in Figure 1 show these precision scores. The thin, colored curves show the true Fisher information, calculated numerically (for method, see Supplementary Appendix H). The thick, black curves show τ̃(x) with v = 2. The two panels in Figure 1 are identical except that Figure 1A has a linear vertical axis and Figure 1B has a logarithmic one. The linear vertical axis facilitates comparison with Clatworthy et al.'s figure 5A, which used linear vertical axes, and the logarithmic vertical axis gives a clearer picture of the data for low r_max values, which are flattened out on the linear axis. The logarithmic vertical axis tends to exaggerate the deviations of the decoding precision from Fisher information because the precision scale is expanded for the worst-matching conditions (those with low precision) relative to the best-matching conditions. 

Figure 1

 

Comparison of contrast decoding precision against the analytical prediction for single Tolhurst-spiking neurons. (A) Squares show Clatworthy et al.'s (2003) Bayesian decoding precision scores for a Tolhurst-spiking neuron with q = 2, c_1/2 = 0.1 (i.e., z = −1), r₀ = 0, and five different r_max values, indicated by different colors; the data were read off from Clatworthy et al.'s figure 5A. Circles show the results of our replication of their simulations. The thin, colored curves plot each neuron's true Fisher information for decoding contrast, calculated numerically (see Supplementary Appendix H). The thick, black curves show the analytical approximation of the decoding precision, τ̃(x), calculated using Equation 13 with K = 1 and v = 2. The vertical black line passing through the peaks of the black curves is positioned along the log contrast axis at log₁₀(c_1/3)—that is, the log of the Michelson contrast for which the mean response is r_max/3. (B) The same as panel A except that the vertical axis has a log scale.

First, note that τ̃(x) is very close to the true Fisher information at the peaks (compare the thick, black lines against the thin, colored lines). This is important because the five observations that we analyze in this section are about the peaks, and our justification for using τ̃(x) to predict decoding precision is that it is approximately equal to the Fisher information. 

For r_max of 50 spikes or more, the horizontal and vertical position of the peak of precision coincides closely with the peak of τ̃(x) (compare symbols against solid curves). This allows us to explain many of Clatworthy et al.'s observations about the precision peak by deriving equations that explain corresponding findings for τ̃(x). For r_max values substantially lower than 50 spikes, the decoding precision does not match τ̃(x) at all well, even at the peak, but Clatworthy et al.'s observations do not apply here either. In the following five subsections, we explain each of the five observations listed previously using the approximation of the decoding precision, τ̃(x), given by Equation 13. 

Observation 1: Multiplicative effect of increasing r_max

From Equation 13, we see that τ̃(x) is proportional to r_max, so increasing r_max increases the precision by the same multiplicative factor for each contrast level; the position of the peak across contrast is unchanged, and the largest absolute change is for the peak precision. This multiplicative change gives rise to the vertical shifts seen for the black curves on the logarithmic vertical axis in Figure 1B. 

Observation 2: Precision peaks slightly below the semisaturation contrast

The exact position of the peak of τ̃(x) along the contrast axis can be found by differentiating τ̃(x) with respect to x, setting the derivative to zero, and solving for x. From Equation 13, we have  

 

Setting this to zero gives  

 

 

 

 

 

Using Equation 15 to substitute for x in Equation 2 gives (for r₀ = 0)  

 

So, regardless of the values of v, r_max, q, or c_1/2, τ̃(x) for a single neuron will always peak at the contrast for which the mean response is r_max/3. We introduce the term c_1/3 for the Michelson contrast that gives rise to a mean response of r_max/3, to be consistent with the term c_1/2 for the semisaturation contrast. The value of log₁₀(c_1/3) is indicated by the black vertical line in each panel of Figure 1, and it passes through the peak of each thick, black curve.  

Observation 3: The peak value of precision is independent of the semisaturation contrast

The peak value of τ̃(x) occurs at a log contrast given by Equation 15. Using this equation to substitute for x in Equation 13, we find that the peak value of τ̃(x) is given by  

 

which is independent of c_1/2.  

Observation 4: The relationship between the maximum precision and q is a steep straight line on log-log coordinates

From Equation 18, we see that max[τ̃(x)] is proportional to q², giving a straight line on log-log coordinates. Figure 2 plots max[τ̃(x)] as a function of q for a single Tolhurst-spiking neuron with r_max = 50 spikes and shows that max[τ̃(x)] is very close to the true peak of precision for this model neuron obtained using Clatworthy et al.'s Monte Carlo methods. 

Figure 2

 

The relationship between the maximum precision and the Naka-Rushton exponent for r_max = 50 spikes. Squares show maximum precision as a function of exponent, q, when r_max = 50 spikes, read off from Clatworthy et al.'s (2003) figure 4C. Circles show the results of our replication of their methods; the differences are negligible. The solid line plots max[τ̃(x)], calculated using Equation 18 with K = 1 and v = 2.

Figure 3

 

Maximum precision as a function of r_max × K. The crosses plot data from Clatworthy et al.'s (2003) figure 5C, showing the height of the peak in decoding precision as a function of r_max × K for a variety of neuronal populations, all with q = 2 and c_1/2 = 0.1. The straight line plots max[τ̃(x)], calculated according to Equation 18 with v = 2.

Observation 5: Linear tradeoff between r_max and number of neurons

Clatworthy et al. (2003) and Chirimuuta and Tolhurst (2005a) both noted that r_max and the number of neurons, K, trade off so that, regardless of the individual values of r_max and K, the decoding precision is a function of their product. In fact, unless r_max × K is very low, the decoding accuracy is close to being proportional to r_max × K. This observation is easily explained by Equation 13, which shows that τ̃(x) is proportional to r_max × K. 

The crosses in Figure 3 replot data from Clatworthy et al.'s (2003) figure 5C. These data show the maximum decoding precision of various populations of neurons, all with q = 2 and c_1/2 = 0.1. The near proportionality between decoding accuracy and r_max × K is indicated by the fact that most of the points lie on a straight line of gradient 1 on these log-log axes. The straight line in our Figure 3 is a plot of max[τ̃(x)], calculated according to Equation 18. It clearly provides a very good match to the decoding precision data for r_max × K ≥ 50. 

Note that although τ̃(x) gives an exact linear trade-off between r_max and K, the true Fisher information for the Tolhurst process does not. As r_max decreases, τ̃(x) tends to underestimate the true Fisher information. 

A more accurate estimate of the decoding precision is given by τ_Tolhurst(x), derived in Supplementary Appendix E:  

 

where  

 

H_Tolhurst(x) is plotted in Figure 4. As the spike rate of each individual neuron increases, H_Tolhurst(x) for that neuron approaches 1/2, and so τ_Tolhurst(x) approaches τ̃(x) (Equation 11) with v = 2. However, as the spike rate decreases, H_Tolhurst(x) increases, and so τ_Tolhurst(x) exceeds τ̃(x) and better reflects the true decoding precision.  

Figure 4

 

H_Tolhurst(x), as defined in Equation 20, plotted as a function of r(x). This function approaches (1 – 1/e) as r(x) → 0 and approaches 0.5 as r(x) → ∞.

Figure 5

 

The inexact nature of the trade-off between r_max and K for the decoding precision of the Tolhurst process. In each panel, the solid line is the same as that in Figure 3 and shows the maximum value of τ̃(x); for this approximation of the decoding precision, the trade-off is exact. The dashed curve shows the maximum value of τ_Tolhurst(x), calculated numerically from Equation 21. Here, the trade-off is only approximate, and the peak of Fisher information does depend slightly on the individual values of r_max and K rather than being simply a function of their product.

For the Naka-Rushton function with r₀ = 0, we can expand Equation 19 in a similar way to Equation 12:  

 

Unlike τ̃(x), τ_Tolhurst(x) is specific to the Tolhurst spiking process rather than being a general approximation of the decoding precision that applies to several different spiking processes.  

In each panel of Figure 5, the solid line is the same as in Figure 3, while the dashed line plots the maximum value of τ_Tolhurst(x). As r_max decreases, τ_Tolhurst(x) and τ̃(x) start to diverge by a factor that approaches (1 – 1/e)/0.5 = 1.264—that is, the ratio of the maximum to minimum values of H_Tolhurst(x). 

This inexact trade-off between r_max and K applies to the Tolhurst process but not to the Poisson or Goris processes: for these, τ̃(x) is derived from the exact expression for the Fisher information, so the trade-off that we derived for τ̃(x) applies to these processes exactly. 

Failure of Clatworthy et al.'s observations for single neurons at low spike rates

When r_max (or, for populations, r_max × K) is substantially less than 50 spikes, the peak of precision does not coincide closely with that of τ̃(x) or even the true Fisher information, so the explanations of Clatworthy et al.'s findings given above are no longer valid. However, most of these findings do not apply for these low r_max values either. The failure of observations 1 and 2 at low spike rates is clear from Figure 1, the failure of observation 5 at low spike rates is shown in Figure 3, and the failure of observation 4 at low spike rates is demonstrated in Figure 6. 

Figure 6

 

The effect of r_max on the relationship between the maximum accuracy and the Naka-Rushton exponent. Results of an investigation the same as that shown in Figure 2 but with three different r_max values (indicated by different colors) and a wider range of exponents, q. Thick, solid lines plot max[τ̃(x)], calculated using Equation 18 with K = 1 and v = 2; symbols (joined by thin lines) plot the height of the peak of decoding precision from Monte Carlo simulations. At high spike rates, the relationship between maximum precision and q is a straight line on log-log coordinates that is well predicted by max[τ̃(x)]; at low spike rates, the relationship between maximum precision and q becomes a shallow curve.

The median r_max for a 200-ms stimulus is only around 5.7 spikes for V1 neurons (Geisler & Albrecht, 1997), suggesting that, with single-neuron models, the spike count has to be implausibly high for the decoding precision to be well approximated by the Fisher information. However, as shown by Clatworthy et al. (2003), Chirimuuta and Tolhurst (2005a), and our Figure 3, r_max can be approximately traded off against the number of neurons so that what is important is the average total spike count of the population rather than that of the individual neurons. With a population code, it is possible to achieve a high total population spike count while keeping the spike counts for the individual neurons at a plausible level. This makes the Fisher information more relevant to understanding population-coding models than coding schemes based on a single neuron. So far we have considered only populations of identical neurons. We now turn to populations of differently tuned neurons. 

Population coding of contrast

Clatworthy et al. (2003) applied maximum-likelihood decoding to three different populations of 18 model Tolhurst-spiking neurons, each with r_max = 10, r₀ = 0, and q = 2. One population had log₁₀(c_1/2) values uniformly distributed between −3 and 0.1, and the other two had c_1/2 values distributed according to the recorded values in either cat or monkey populations, found by arranging the neurons in ascending order of c_1/2 and then sampling the population at equal percentile intervals. We did not have access to the exact sets of cat or monkey c_1/2 values that they used, but we estimated them by fitting functions to the cat or monkey c_1/2 distributions given in Clatworthy et al.'s figure 6 and then sampling these distributions in equal percentile steps (see Supplementary Appendix I for details). The advantage of our method is that it can easily be extended to neuronal populations of any size. Having set up the populations of neurons, we then calculated decoding precision scores using Clatworthy et al.'s Monte Carlo methods (see Supplementary Appendix G for details). Figure 7 shows that our decoding precision scores are very similar to those of Clatworthy et al., confirming that our method of generating the sets of c_1/2 values is a good approximation to that of Clatworthy et al. The colored curves in Figure 7 show τ̃(x), calculated using Equation 12 with v = 2. Since r_max = 10 in these simulations, H_Tolhurst(x) ≈ 0.5 for the most informative neurons, so τ_Tolhurst(x) ≈ τ̃(x); there was no advantage in using the closer but more complex approximation. Figure 7 shows that, even for these small populations of neurons, with r_max of only 10 spikes, the decoding precision from the Monte Carlo simulations is very close to τ̃(x) over a wide range of contrasts. 

Figure 7

 

Population decoding precision for three different populations of K = 18 neurons. Within each population, the log₁₀(c_1/2) values were uniformly distributed between log contrasts of −3 and 0.1 (red) or distributed according to recorded cat (green) or monkey (blue) populations (see Supplementary Appendix I for details). For each neuron, r_max = 10, r₀ = 0, and q = 2. Squares show Clatworthy et al.'s (2003) decoding precision data, read off from their figure 7. Circles show the decoding precision data from our replication of their Monte Carlo methods. The solid curves show τ̃(x), calculated using Equation 12 with v = 2. Note that because our method of selecting the c_1/2 values differed slightly from that of Clatworthy et al., our decoding precision scores show a very slight systematic difference from theirs, which is more noticeable in the cat data. τ̃(x) was calculated using our c_1/2 values and therefore shows a slightly better match to our precision scores than to those of Clatworthy et al.

In Figure 7, the only substantial deviation of decoding precision from Fisher information occurs at the ends of the contrast range, where the precision scores shoot up. This is an artifact, discussed by Clatworthy et al., caused by the fact that in the simulations the likelihood functions were calculated over a finite range of contrasts. This means that on trials where the inferred contrast would have fallen beyond the ends of the contrast range, the inferred contrasts instead pile up on the ends of the range, artificially boosting the number of trials on which the inferred contrast takes those values. Thus, when the stimulus contrast really is at or close to one of the ends of the range, it is likely to be close to the inferred contrast, so the precision is artificially high. The Fisher information does not show this effect because it is a purely local measure, derived from the contrast-response functions at each point along the contrast axis, so it is unaffected by any bounds on the range of contrasts. 

τ̃(x) is a good approximation of the decoding precision for the Tolhurst process as long as r_max is not too low. For low values of r_max, τ̃(x) tends to underestimate the true Fisher information. As noted earlier, this means that there is not an exact trade-off between r_max and the number of neurons, K. For low values of r_max, it is better to use τ_Tolhurst(x) to accurately predict decoding precision. This is demonstrated in Figure 8. 

Figure 8

 

Results of four sets of simulations the same as those in Figure 7 except with different values of r_max and K, indicated at the top of each panel. As in Figure 7, circles plot decoding precision from the Monte Carlo simulations and solid curves plot τ̃(x), calculated using Equation 12 with v = 2. Dashed curves plot τ_Tolhurst(x), calculated using Equation 21. For r_max = 16, τ̃(x) and τ_Tolhurst(x) are very well matched to each other and to the true decoding precision. For r_max = 4 and below, only τ_Tolhurst(x) is well matched to the true decoding precision. This figure also illustrates the breakdown of the linear trade-off between r_max and K when r_max is low. Each panel has the same product r_max × K, but the decoding precision falls with increasing r_max. This breakdown of the linear trade-off applies to the Tolhurst process shown here but not to the ordinary Poisson process or the Goris process.

There are at least two reasons why we sometimes need an expression for the decoding precision that remains a close approximation down to very low spike rates. First, the median spike count for a 200-ms stimulus presentation is only around 5.7 spikes for a V1 neuron tuned to the stimulus (Geisler & Albrecht, 1997), and observers can make complex judgements based on exposures even shorter than that (Thorpe, Fize, & Marlot, 1996). Second, in many situations the majority of neurons will not be well tuned to the stimulus but may still contribute to task performance due to their abundance; these neurons would be expected to have a very low average spike count over the course of a stimulus presentation. 

The psychometric function for contrast detection

Chirimuuta and Tolhurst (2005a) used their contrast coding model to simulate a two-alternative forced-choice (2AFC) contrast detection task, in which the observer is presented with a zero-contrast nontarget stimulus and an above-zero-contrast target stimulus and has to pick the target. Chirimuuta and Tolhurst (2005a) measured the model's proportion of correct responses as a function of the target contrast and fitted a Weibull psychometric function to the data. 

For 2AFC tasks, the Weibull function can be defined as  

 

α is the threshold—that is, the target contrast that gives P(correct) = 1 – 0.5/e = 0.816…—and β controls the function's shape on linear axes or slope on log axes (see May & Solomon, 2013, for an in-depth analysis of the Weibull function). In psychophysical contrast detection experiments with human observers, β usually takes a value of about 3 (Foley & Legge, 1981; Nachmias, 1981; Mayer & Tyler, 1986; Meese, Georgeson, & Baker, 2006; Wallis, Baker, Meese, & Georgeson, 2013).  

When Chirimuuta and Tolhurst's (2005a) model had a standard Naka-Rushton contrast-response function with r₀ = 0 and q = 2, Weibull β for detection varied from 1.75 to 1.99 as the number of neurons increased from 1 to 23, but the model's β never reached the normal human level of around 3. Chirimuuta and Tolhurst then introduced a threshold to the Naka-Rushton function by subtracting a constant value from the output and setting negative values to zero. With a threshold on the Naka-Rushton function, β ranged from 2.25 to 4.20, providing a better match to psychophysically obtained values. Chirimuuta and Tolhurst assumed that their failure to obtain sufficiently high Weibull β values with the standard Naka-Rushton function had been caused by the lack of a threshold, and they suggested that the standard, unthresholded Naka-Rushton function “may be crucially inadequate” as a model of the neuronal contrast-response function (Chirimuuta & Tolhurst, 2005a, p. 2956). However, we now show that, when r₀ = 0, the model's psychometric function is close to a Weibull function with β = q. Thus, the real reason that Chirimuuta and Tolhurst always obtained a β of about 2 with the standard Naka-Rushton function is that they always kept q at 2 in these simulations. We show that one can obtain any Weibull β by setting q close to the required β value. 

Proof that, when r₀ = 0, the model's psychometric function for 2AFC contrast detection is close to a Weibull function with β = q

At the low contrasts corresponding to the model's contrast detection threshold, the population spike rate is typically so low that the Fisher information does not provide a useful approximation of performance. Therefore, we need to use other methods. For contrast detection, the model's performance can be derived straightforwardly from basic probability theory. 

Because r₀ = 0, there is zero response to zero contrast: The nontarget stimulus can never elicit a single spike. If the target elicits at least one spike, the model will respond correctly. If the target fails to elicit any spikes, the model has to guess and will be correct half the time on a 2AFC task. In summary, the model will be correct on all 2AFC trials except half of those on which the target failed to elicit any spikes. This statement can be formalized as follows:  

 

where P(no spikes) is the probability of getting no spikes in response to the target. We can already see that the model's psychometric function has a similar form to the Weibull function: Display Formula  

From Equation E.4 of Supplementary Appendix E, we see that, for a single Tolhurst-spiking neuron,  

 

Note that since the Weibull function is defined as a function of contrast in Michelson (linear) units, we express the mean response, r(·), as a function of Michelson contrast (rather than log contrast) in Equation 24. For a population of K statistically independent neurons,  

 

where r_j(c) is the contrast-response function of neuron j. Using Equation 1 to expand r_j(c), we have  

 

where (r_max)_j and (c_1/2)_j are the r_max and c_1/2 parameters, respectively, of neuron j. Using Equation 26 to substitute for P(no spikes) in Equation 23, we get an exactly correct expression for the model's psychometric function for contrast detection:  

 

If all the neurons in the population being monitored by the observer have the same contrast-response function, then Equation 27 reduces to  

 

In this case, Equation 28 shows that there is an exact linear trade-off between r_max and K in the psychometric function for 2AFC detection rather than the approximate trade-off that we get with the Fisher information.  

Equations 27 and 28 are precisely correct, but it is illuminating to derive an approximation. As long as there are enough neurons or r_max is sufficiently high, the contrast detection threshold will be somewhat lower than the lowest c_1/2 in the population. Thus, at threshold (i.e., around the middle of the psychometric function), the denominators of Equations 26 to 28 become dominated by the Display Formula term, and the c^q term will make little difference. Dropping c^q from the denominator of Equation 27, we get  

 

where α is a constant, given by  

 

Relation 29 has the form of a Weibull function with β = q. The near equality in Relation 29 approaches equality as K or r_max increase. □  

The derivation of Relation 29 is similar to the derivation of the psychometric function for Quick's (1974) vector-magnitude model of contrast detection. The performance of the vector-magnitude model is exactly equivalent to that of a model in which the observer monitors a number of detectors that each have a detection probability that can be described by a Weibull function, with the same β for all detectors. Nachmias (1981) called this assumption of identical β for each detector the “homogeneity assumption,” and it is largely equivalent to the implicit assumption in the above proof that all neurons have the same q. In the Discussion we expand on the links between our analysis and that of Quick (1974). 

Another attribute of the model's psychometric function: Lapse rate

Note that Relation 29 is a Weibull function with zero lapse rate; that is, it predicts that the model's performance will asymptote to perfect performance, P(correct) = 1, with increasing target contrast. For most instantiations of the model, this is very close to the truth. However, when r_max and K are both very low, the model's asymptotic performance is far below 1. This is because as c increases, the denominator of Equation 26 becomes more and more dominated by the c^q term, and the c_1/2^q term makes less and less difference, so  

 

For very low r_max and K, this value can be substantially above zero so that, even for infinite contrast, the model has to guess on a significant proportion of 2AFC trials.  

This behavior can be accommodated by including a “lapse rate” parameter, λ, in the definition of the Weibull function:  

 

 

This equation approaches an asymptote of P(correct) = ( 1 – λ) as c → ∞. When λ = 0, Equation 32 reduces to Equation 22. We now derive an expression for the model's “lapse rate” parameter. (Strictly speaking, the model never lapses; a low asymptotic performance level arises from a low spike rate rather than a finger error or a failure to look at the stimuli on some 2AFC trials.) 

Using Equation 31 to substitute for P(no spikes) in Equation 23, we obtain the asymptotic value of P(correct) given by  

 

Since this asymptotic value of P(correct) is (1 – λ), we have  

 

This expression for λ quickly approaches zero as r_max or K increase.  

Verifying our equations using Monte Carlo simulations

We simulated the 2AFC contrast detection experiments for a range of model parameterizations (see Supplementary Appendix J for details of the method). For each parameterization, each neuron in the population had an identical contrast-response function. We fitted the three-parameter Weibull function (Equation 32) to the data for each model parameterization. The fitted values of α, β, and λ are plotted as symbols in Figures 9, 10, and 11, respectively. The solid lines plot the corresponding analytical expressions (α given by Equation 30, β given by q, and λ given by Equation 34). 

Figure 9

 

Filled circles show Weibull α parameters fitted to the 2AFC contrast detection simulation data. Each column of panels shows a different r_max; each row shows a different q. The solid lines plot the approximation of α given by Equation 30.

Figure 10

 

Symbols show Weibull β parameters fitted to the 2AFC contrast detection simulation data. Each panel shows a different r_max. Different values of q are indicated by different symbols and colors according to the legend on the right. The solid lines plot the predicted values of β (= q) for our approximation of the psychometric function (Relation 29).

Figure 11

 

Filled circles show Weibull λ parameters fitted to the 2AFC contrast detection simulation data. Each column of panels shows a different r_max; each row shows a different q. The solid curves plot λ given by Equation 34. As shown in Equation 34 (and borne out by our simulations), λ depends only on K and r_max and is independent of q.

Equation 34 is the model's true lapse rate parameter (not an approximation), so it is not surprising that it matches the fitted λ values extremely well in Figure 11. Our analytical expressions for α and β are approximations that become increasingly accurate as K or r_max increase. 

The Consul-Jain spiking process also generates a Weibull psychometric function for detection

The analytical expressions for the psychometric function for detection, derived above, apply only to the Tolhurst spiking process. We can also derive analogous expressions for the Consul-Jain process (which includes the ordinary Poisson, when F = 1). Equation E.17 of Supplementary Appendix E states that, for a single Consul-Jain-spiking neuron,  

 

where F is the Fano factor. If we follow a series of mathematical steps analogous to Equations 23 to 29 above but use Equation 35 instead of Equation 24 to express the probability of a single neuron not spiking, we obtain an approximation of the psychometric function with the same form as Relation 29 but with  

 

where F_j is the Fano factor or neuron j. Similarly, if we follow a series of steps analogous to those in Equations 31 to 34 but use Equation 35 instead of Equation 24 to express the probability of a single neuron not spiking, we obtain the following expression for the “lapse rate” parameter:  

 

 

Discussion

The purpose of this study was to explain a number of empirical modeling results reported by Tolhurst and colleagues. Their results were obtained from Monte Carlo simulations using models of spiking neurons with Naka-Rushton contrast-response functions and the doubly stochastic Poisson spiking process defined in Equation 5. The numerical nature of the simulations meant that it was not clear why the results occurred or how they generalized across parameter space. 

We addressed these problems by deriving equations to explain the model's performance. This kind of analysis is essential if we are to understand the brain: If we use realistic models to simulate brain processes but do not understand why the models behave in the way that they do, then we have not really explained anything; we have just shifted the problem from understanding the brain to understanding the model. 

We began by deriving a closed-form expression for the Tolhurst likelihood function that was more tractable than the infinite series that originally defined this function. This expression played a role in understanding two facets of the model's performance: (a) the relationships between decoding precision and the neuronal parameters, and (b) the form of the psychometric function for contrast detection. 

Decoding precision

To explain how the neuronal parameters map onto decoding precision, we derived an analytical approximation of the decoding precision, τ̃(x), which can be adjusted to apply to a variety of different spiking processes by setting the value of a single scalar parameter, v. For the Tolhurst and Consul-Jain processes, v should be set to the Fano factor. For the Goris process, v should be set to 1/1(1 − Display Formula ), where Display Formula is the variance of the gain signal. For the Tolhurst and Consul-Jain processes, τ̃(x) is an estimate of the Fisher information. For the Goris process, the Fisher information varies across trials, and τ̃(x) is its modal value.  

The expression for τ̃(x) revealed some surprisingly simple relationships between decoding precision and the neuronal parameters, and explained the five observations of Clatworthy et al. (2003) that we investigated. For example, Equation 18 shows that, for a population of identical, statistically independent neurons, the height of the peak of decoding precision is approximately proportional to r_max × K × q² and is independent of c_1/2. The expression for τ̃(x) also revealed that, to a good approximation, the contrast most accurately encoded by the neuron is that for which the mean response is r_max/3; we call this contrast c_1/3. Figure 7 shows that τ̃(x) matches the decoding precision very closely for a population of only 18 Tolhurst-spiking neurons with r_max of only 10 spikes. 

τ̃(x) is an exact expression for the modal Fisher information for the Goris process and gives the exact expression for the Fisher information for the Poisson process, which is the Consul-Jain process with a Fano factor of 1. However, for the Tolhurst process and the Consul-Jain process with F > 1, τ̃(x) is an approximation of the Fisher information and is less accurate when the mean spike rate of the most informative neurons is very low. In Supplementary Appendix E, we derived expressions, τ_Tolhurst(x) and τ_C-J (x), that are close approximations of the Fisher information of the Tolhurst and Consul-Jain processes across all parameter values. These expressions reveal more complicated relationships between decoding precision and the neuronal parameters that hold when the spike rate is very low. Figure 8 shows the superiority of τ_Tolhurst(x) over τ̃(x) at very low spike rates. 

The psychometric function for contrast detection

To explain Chirimuuta and Tolhurst's (2005a) finding regarding the slope of the psychometric function for contrast detection, we derived an analytical approximation of the model's psychometric function, which we showed approaches a Weibull function with β = q as r_max or K increase. We thus refuted Chirimuuta and Tolhurst's conclusion that the standard Naka-Rushton function is unable to give Weibull β values that are high enough to match those of human observers. To obtain β values of around 3, typical of human observers, we need a Naka-Rushton exponent of around 3. Such levels are not atypical, and indeed Table 1 shows that the mean q over 628 cells from three physiological studies of V1 is 2.9. 

Table 1

 

Mean and median Naka-Rushton exponent, q, from three physiological studies of V1. Albrecht and Hamilton (1982) did not report the median, so these are omitted. Sclar, Maunsell, and Lennie (1990) did not report the mean, so the value given in this table was calculated from the histogram of q values shown in their figure 6. The centers of the bars of their histogram ranged from 0.875 to 6.375 in steps of 0.25. The numbers of neurons in the corresponding bins were 8, 5, 2, 12, 7, 6, 4, 3, 9, 3, 3, 5, 4, 4, 0, 4, 1, 1, 1, 2, 0, 0, and 1. Taking the q of each neuron in a histogram bin to be the central value of the bin (a reasonable approximation given the narrowness of their bins), this gives a mean value of 2.65, calculated from 85 neurons. Sclar et al.'s (1990) study included one further V1 neuron, with q > 8, but because they were no more specific about its q value than that, we excluded it from our analysis. The overall mean was the average of the mean q values, weighted by the number of neurons in each study (this is equivalent to pooling the neurons across all three studies and finding the mean of the pooled sample).

 

This example illustrates the power of the analytical approach taken here, compared with Monte Carlo simulations. Because of the considerable time that it takes to run the simulations, Chirimuuta and Tolhurst (2005a) could not feasibly explore every corner of parameter space in each investigation. Therefore, they took the entirely reasonable step of fixing the Naka-Rushton exponent, q, at a physiologically plausible level of 2. However, this turned out to be a fatal decision because it prevented them from ever finding a parameterization of the model that gave a sufficiently high Weibull β with the standard Naka-Rushton tuning function. Using the analytical approach, Relation 29 makes it immediately clear that the Naka-Rushton exponent, q, is the key parameter for controlling β, and that the model's psychometric function approximates a Weibull function with β actually equal to q. 

It should be noted that our derivation of a psychometric function with the form of a Weibull function applies only to the case of r₀ = 0. With nonzero r₀, the analytical form of the psychometric function is different (we have analyzed this more general case and will present it in another article). It is implausible that, in human vision, contrast detection is mediated entirely by neurons with zero spontaneous firing rate. The assumption that there is no neuronal response to zero contrast is often called the high-threshold assumption.² Under the conventional assumption of additive, stimulus-independent noise, the lack of response to zero contrast implies that there is a threshold on the output of the sensory units that lies enough standard deviations above the mean of the noise for there to be a negligible probability of a response to zero contrast. Because Chirimuuta and Tolhurst's model has no sensory response to zero contrast, it is formally equivalent to a high-threshold model even though (with the standard Naka-Rushton function) it does not actually contain a sensory threshold. In high-threshold theory (and in Chirimuuta and Tolhurst's model), detection errors are always unlucky guesses on 2AFC trials that failed to elicit a response, whereas there is plenty of psychophysical evidence that incorrect responses are caused at least partly by hallucinations due to noise in the nontarget interval rather than entirely by unlucky guesses (Tanner & Swets, 1954; Swets, Tanner, & Birdsall, 1961; Nachmias, 1981; Solomon, 2007; Laming, 2013). 

Despite the implausibility of Chirimuuta and Tolhurst's model with r₀ = 0 as a model of contrast detection, we presented our analysis of it for three reasons. 

1.  
   It allows us to fully understand why Chirimuuta and Tolhurst's (2005a) contrast detection simulations resulted in a fitted Weibull β that approached 2 with increasing number of neurons. This in turn allows us to refute their conclusion that the Naka-Rushton function requires a threshold to make it a plausible model of the neuronal contrast-response function.
2.  
   Neurons with zero or negligible spontaneous firing rates do exist (e.g., see Dean, 1981a, figure 2), so it is not inconceivable that there are some organisms or experimental situations to which our analysis applies.
3.  
   It allows us to address Tyler and Chen's (2000) claim that high-threshold analysis of probability summation is “fundamentally flawed.” This idea is explored in the next subsection.

High-threshold probability summation

Our derivation of the model's psychometric function for detection is an example of high-threshold probability summation. When r₀ = 0, each neuron acts as a detector; the observer detects the target if at least one neuron responds during the target presentation, and has to guess the correct answer otherwise. The more neurons the observer is monitoring, the greater the chance that at least one neuron will respond. The term probability summation refers to this increase in detection probability due to an increase in the number of detectors. The psychometric function in this case gives the probability that at least one neuron responds during the target presentation, or the observer guesses correctly if none respond. 

Quick (1974) showed that if the detectors are statistically independent and each detector has a detection probability that is a Weibull function with the same β, then the observer's psychometric function will be a Weibull function with that β-value and with the detection threshold parameter, α, given by  

 

where α_j is the detection threshold parameter of detector j. For a very clear derivation of Equation 38, see Nachmias (1981), but note that his equation 4 has a typographical error: It is missing the minus sign on the exponent, β. If all the detectors are identical (but statistically independent), then they all have the same α_j, and Equation 38 reduces to  

 

 

Thus, in Quick's model, the Weibull β parameter controls how much of a reduction in detection threshold we achieve by increasing the number of detectors, K. If β is low, then the detection threshold, α, decreases quickly as K increases; if β is high, then the detection threshold decreases more slowly as K increases. 

Equation 38 gives the observer's detection threshold according to Quick's model, while Equation 30 gives the observer's approximate detection threshold according to Chirimuuta and Tolhurst's (2005a) model. Equation 38 has the exact form of Equation 30 if β = q and  

 

Thus, there is a near equivalence between Quick's model and that of Chirimuuta and Tolhurst. If all the neurons in Chirimuuta and Tolhurst's model have the same contrast-response function, then the model's detection threshold approximation given by Equation 30 reduces to  

 

Therefore, since β ≈ q in Chirimuuta and Tolhurst's model, their model shows approximately the same probability summation effects as Quick's model, with detection threshold proportional to K^−1/β. For the modeling in Figure 9, all the neurons were identical, so Equation 41 is equivalent to Equation 30 (which was used to generate the solid lines in Figure 9), and it is clear that this equation does accurately predict the detection threshold of Chirimuuta and Tolhurst's model, particularly for the higher values of r_max or K.  

As noted previously, there is psychophysical evidence against Quick's high-threshold model. Tyler and Chen (2000) went further, arguing not just that high-threshold theory has evidence against it but that the analysis of probability summation using high-threshold theory is “fundamentally flawed.” Their analysis of high-threshold theory made several basic assumptions about the high-threshold model, including the following five. 

1.  
   The observer monitors a set of channels.
2.  
   Within each channel is a continuous signal that increases linearly with the stimulus strength.
3.  
   Noise is added to this signal.
4.  
   The noise might be additive (standard deviation independent of the signal level) or multiplicative (standard deviation proportional to a power function of the signal level), or the noise could be the sum of an additive and a multiplicative component.
5.  
   A sensory threshold is applied to the noisy internal signal in each channel so that a stimulus is detected if and only if the noisy signal falls above threshold in at least one of the channels.

We should emphasize that Tyler and Chen (2000) were not arguing in support of this model; they argued that it was fundamentally flawed. 

Tyler and Chen noted that if the internal noise in this model were fully multiplicative, then the internal signal-to-noise ratio would be unchanged by a change in stimulus strength. Thus, detection performance would remain the same for any stimulus level, and measurement of detection thresholds would be impossible! They concluded that even if there is a multiplicative component to the noise, detection performance must be limited by an additive component. 

Tyler and Chen derived the shape of the additive noise distribution that would give rise to a Weibull psychometric function for the high-threshold model outlined previously.³ They showed that, for most Weibull β values, the noise distribution deviated markedly from a Gaussian. They found this unacceptable because, according to the central limit theorem, the sum of a large number of non-Gaussian-distributed random variables is asymptotically Gaussian distributed. Therefore, if there are many different sources of noise from the stimulus to the neuronal decision mechanism, the noise is likely to be Gaussian at the decision mechanism. They then showed that high-threshold probability summation fails for additive Gaussian noise. They showed that, if the sensory threshold were low enough to be exceeded by the noisy signal in one channel 75% of the time, then the noisy signal would exceed the sensory threshold in at least one of 100 channels almost all the time, even when the stimulus intensity was reduced to zero. Thus, if the stimulus area or number of components increased so that the observer was monitoring many more channels, the observer would almost always be in a detect state, even when the stimulus was absent. The signal would have to be reduced to a physically unachievable negative contrast for the observer to be in a detect state 75% of the time. This problem does not occur with the noise distribution implied by the Weibull function because the Weibull probability density function falls to zero at the sensory threshold (see Supplementary Appendix K, especially Figure K.1); therefore, when the stimulus intensity is zero, no detector's noisy signal will exceed its sensory threshold. However, as already noted, Tyler and Chen (2000) rejected that distribution because of its “bizarre” non-Gaussian nature when β is not close to 4 (p. 3127). They therefore concluded that probability summation with high-threshold theory is fundamentally flawed. 

Tyler and Chen's arguments are perfectly valid, but they apply only to the set of high-threshold models defined by their assumptions (i.e., that each channel contains a continuous signal to which noise is added). With a detection model based on spiking neurons, the signal is discrete. When the spontaneous firing rate is zero, as it was in Chirimuuta and Tolhurst's model, detection occurs when a single neuron produces a single spike. In these circumstances, it is not appropriate to approximate the neuronal signal as a continuous signal to which noise is added. It is perfectly legitimate to apply probability summation to find the probability that at least one neuron gives at least one spike. When we do this, we find that the model's threshold is very close to being proportional to K^−1/β, as in standard high-threshold Weibull analysis. The contrast threshold for detection never reaches zero regardless of how many neurons the observer is monitoring, but we do not have to rely on bizarre model characteristics to achieve this. Aside from the zero spontaneous firing rate, the neurons in Chirimuuta and Tolhurst's model have contrast-response functions and noise distributions that are physiologically plausible to a reasonable extent. One could ask what the equivalent “continuous signal plus additive noise” model is. Since Chirimuuta and Tolhurst's model produces a psychometric function that closely approximates a Weibull function, the equivalent “continuous signal plus additive noise” model is closely approximated by the one derived by Tyler and Chen (and outlined in Supplementary Appendix K), with the “bizarre” noise distributions. However, this bizarreness comes from forcing Chirimuuta and Tolhurst's model into the Procrustean bed of “continuous signal plus additive noise” rather than being an implausible characteristic of the model itself. 

Conclusions

We derived equations that explained the performance of the contrast coding model described by Tolhurst and colleagues (Clatworthy et al., 2003; Chirimuuta & Tolhurst, 2005a). These equations gave a deep insight into the results of their Monte Carlo simulations. 

As long as the stimulus contrast is high enough (and the neuronal population is large enough) to generate a sufficiently high population spike rate, the decoding precision can be closely approximated by the Fisher information. We derived an estimate of decoding precision, τ̃(x), which approximates the Fisher information for a population of neurons with Tolhurst's spiking process as long as the mean spike count of the most informative neurons is around five spikes or more. τ̃(x) is also an estimate of the Fisher information for the Consul-Jain spiking process. Furthermore, it gives the exact Fisher information for the Poisson process and the exact modal value of the Fisher information for the Goris process when the decoder has access to the gain signal in Goris et al.'s (2014) model of neuronal spiking. Our expression for τ̃(x) revealed simple relationships between the properties of the neurons and the decoding precision that hold to a good approximation when the mean count rate is sufficiently high. We used this expression to explain five relationships between decoding precision and the neuronal parameter values that Clatworthy et al. (2003) observed from their Monte Carlo simulations. 

For the case of contrast detection, the spike rate is too low for the Fisher information to match decoding precision. To analyze the performance of Chirimuuta and Tolhurst's (2005a) model in a 2AFC contrast detection task, we used basic probability theory. We derived an expression for the model's psychometric function for contrast detection and showed that as K or r_max increase, the psychometric function asymptotically approaches a Weibull function with β = q. Our work therefore reveals a previously unknown theoretical connection between two of the most widely used functions in vision science: the Weibull psychometric function and the Naka-Rushton contrast-response function. This relationship explained why Chirimuuta and Tolhurst always obtained a Weibull β of about 2 in their modeling (they always had q = 2 in their assessments of the model's Weibull β) and allowed us to refute their conclusion that it is necessary to have a threshold on the Naka-Rushton function to achieve Weibull β values that match those found with human observers. Their threshold on the Naka-Rushton function had a similar effect to increasing q, as it made the spike rate increase more abruptly with increasing contrast. 

Acknowledgments

This work was supported by EPSRC grant EP/H033955/1 to Joshua Solomon. The authors thank Christopher Tyler for giving detailed feedback on an earlier draft of this article. We also thank Wei Ji Ma and Peter Neri for insightful reviews of our manuscript. 

Commercial relationships: none. 

Corresponding author: Keith A. May. 

Email: [email protected]. 

Current address: UCL Department of Computer Science, London, UK. 

Footnotes

¹  In this paper, we use the word trial in two ways. First, we use it in the way a physiologist would, to mean a stimulus presentation. Second, we use it to mean a trial on a two-alternative forced-choice (2AFC) psychophysical experiment, in which the observer is presented with two stimuli and has to make a response. To distinguish these two meanings, we always refer to the latter type of trial as a 2AFC trial.

Footnotes

²  In this discussion, we use the word threshold in two ways: (a) to refer to an internal threshold on the sensory signal, and (b) to refer to the stimulus contrast corresponding to a particular level of detection performance. We have tried to make the meaning clear in each case by using the term sensory threshold for the former case and the term detection threshold for the latter.

Footnotes

³  As pointed out by Mortensen (2002), Tyler and Chen's (2000) published equation for the probability density function (PDF) of the noise (Tyler and Chen's equation 4b) contains errors. However, Tyler and Chen's plots of the noise PDFs (shown in their figure 2b) are correct, so we assume that the errors in Tyler and Chen's equation 4b are typographical errors rather than fundamental problems with their analysis. To clarify matters, we present in Supplementary Appendix K a derivation of the PDF that is based on Mortensen's derivation but is hopefully easier to follow than Mortensen's derivation or that of Tyler and Chen.