In a two-choice psychophysical task, such as a stimulus discrimination or detection task, a neuron is said to contain a ‘choice-related signal’, or ‘decision-related signal’ when its activity carries information about the behavioral choice above and beyond the information that it carries about the stimulus (Britten et al., 1996; Parker and Newsome, 1998; Nienborg et al., 2012). The interpretation of choice-related signals in terms of decision-making mechanisms is however difficult. Much progress in our understanding of their meaning has relied on using models to derive mathematically the relationship between the underlying decision-making mechanisms and different measures of activity-choice covariation (Haefner et al., 2013; Pitkow et al., 2015) usually used to quantify choice-related signals.

The most widely used measure of activity-choice covariation for tasks involving two choices is choice probability, $\mathrm{CP}$. The $\mathrm{CP}$ is defined as the probability that a random sample of neural activity from all trials with behavioral choice $D$ equal to 1 is larger than one sample randomly drawn from all trials with choice $D=-1$ (Britten et al., 1996; Parker and Newsome, 1998; Nienborg et al., 2012; Haefner et al., 2013):

where $r$ is any measure of the neural activity, which we will here consider to be the neuron’s per-trial spike count. Another prominent measure of choice-related signals is choice correlation, $\mathrm{CC}$ (Pitkow et al., 2015). This quantity is defined under the assumption that the binary choice $D$ is mediated by an intermediate continuous decision value, $d$. This value may represent the brain’s estimate of the stimulus, or an internal belief about the correct choice. The definition of CC further assumes that the categorical choice $D$ is related to $d$ via a thresholding operation such that the choice depends on whether $d$ is smaller or larger than a threshold θ (Gold and Shadlen, 2007). Its expression is as follows:

where $\mathrm{cov}(r,d)$ is the covariance of the neural responses with $d$, and $\mathrm{var}r$, $\mathrm{var}d$ their variance across trials. Perhaps, the simplest measure of activity-choice covariation, which has been used in empirical studies (Mante et al., 2013; Ruff et al., 2018), is what we called the choice-triggered average, $\mathrm{CTA}$, defined as the difference between a neuron’s average spike count $r$ across trials with behavioral decision $D=1$ minus the average spike count in trials with decision $D=-1$:

The CP and CTA quantify activity-choice covariations without assumptions about the underlying decision-making mechanisms. However, their interpretation has commonly (Nienborg et al., 2012) been informed in previous analytical and computational studies by assuming a specific feedforward decision-threshold model of choice-related signals (Shadlen et al., 1996; Cohen and Newsome, 2009b). Haefner et al., 2013 used that model to derive an analytical expression for CP valid under two assumptions that are often violated in practice: first, the model assumes a causally feedforward structure in which sensory responses caused the decision, and second, it is assumed that both decisions are equally likely. However, the presence of informative stimuli leads to one choice being more likely than the other, hampering the application of the analytical results to Grand CPs and to detection tasks (Bosking and Maunsell, 2011; Smolyanskaya et al., 2015), which involve informative stimuli. Furthermore, decision-related signals have empirically been shown to reflect substantial feedback components (Nienborg and Cumming, 2009; Nienborg et al., 2012; Macke and Nienborg, 2019). We will next extend this previous model (Haefner et al., 2013) to obtain a general expression of the CP valid for informative stimuli and regardless of the feedforward or feedback origin of the dependencies between the neural responses and the decision variable.

We first consider a most generic model in which we simply assume that the response r_i of the $i-th$ sensory neurons covaries with the behavioral decision $D$, but without making any assumption about the origin of that covariation (Figure 1A). We find that to a first approximation (exact solution provided in Methods), the CP of cell $i$ captures the difference between the distributions $p(r_i|D=1)$ and $p(r_i|D=-1)$ resulting from a difference in their means, and hence is related to the CTA:

Figure 1

Download asset Open asset

The $\mathrm{CTA}$ generically quantifies the linear dependencies between responses and choice, and this approximation of the CP does not depend on their feedforward or feedback origin (Figure 1A). We next add the assumption that the relationship between a neuron’s response and the choice is mediated by the continuous variable $d$, as commonly assumed by previous studies and described above (Figure 1B). This splits any correlation between the neural response r_i and choice $D$ into the product of the two respective correlations: ${\displaystyle \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}(r_i,D)=\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}(r_i,d)\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}(d,D)}$ = ${\displaystyle {\mathrm{C}\mathrm{C}}_i\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}(d,D)}$, where ${\mathrm{CC}}_i=\mathrm{corr}(r_i,d)$ is the choice correlation as defined in Equation 2. It follows (see Methods) that:

where ${\mathrm{CTA}}_d$ is the average difference in $d$ between the two choices, in analogy to the ${\mathrm{CTA}}_i$ for neuron $i$. Equation 5 describes how activity-choice covariations appear in the model (Figure 1C): the threshold mechanism dichotomizes the space of the decision variable, resulting in a different mean of $d$ for each choice, which is quantified in ${\mathrm{CTA}}_d$. If the activity of cell $i$ is correlated with the decision variable $d$ (non zero ${\mathrm{CC}}_i$), the ${\mathrm{CTA}}_d$ is then reflected in the ${\mathrm{CTA}}_i$ of the cell. In previous theoretical work (Haefner et al., 2013), the distribution over $d$ was assumed to be fixed and centered on the threshold value θ. Here, we remove that assumption and consider that $d$ may not be centered on the threshold if the stimulus is informative, containing evidence in favor of one of the two choices, or if the choice is otherwise biased. In those cases, the normalized ${\mathrm{CTA}}_d$ in Equation 5, namely ${\mathrm{CTA}}_d/\sqrt{\mathrm{var}d}$, can be determined (see Materials and methods) in terms of the probability of choosing choice 1, $p_{\text{CR}}\equiv p(D=1)=p(d>\theta )$, which we call the ‘choice rate’, $p_{\text{CR}}$. Since the decision variable is determined as the combination of the responses of many cells, its distribution is well approximated by a Gaussian distribution, but now with a nonzero mean determined by the stimulus content. With this assumption, the normalized ${\mathrm{CTA}}_d$ for $p_{\text{CR}}=0.5$ is equal to $4/\sqrt{2\pi }$, and for each other $p_{\text{CR}}$ value differs by a scaling factor

where $\varphi (x)$ is the density function of a zero-mean, unit variance, Gaussian distribution, and ${\mathrm{\Phi }}^{-1}$ is the corresponding inverse cumulative density function. By construction, $h(p_{\text{CR}})=1$ for $p_{\text{CR}}=0.5$ where it has its minimum. Given the factor $h(p_{\text{CR}})$, combining Equations 4 and 5 we can relate CP and CC across different ratios $p_{\text{CR}}$, corresponding to different stimulus levels, irrespectively of whether CP is caused by feedforward or feedback signals. In the linear approximation (see Methods for the exact formula and derivation with the decision-threshold model), this relationship reads:

For equal fractions of choices, $p_{\text{CR}}=0.5$, this CP expression corresponds to the linear approximation derived in Haefner et al., 2013. Note that extending the CP formula to $p_{\text{CR}}\ne 0.5$ required us to also make explicit the dependency of the choice correlations on the choice rate, ${\displaystyle {\mathrm{C}\mathrm{C}}_i(p_{\mathrm{C}\mathrm{R}})}$. Unlike $h(p_{\text{CR}})$ which is an effect of the decision-making threshold mechanism and shared by all neurons, ${\displaystyle {\mathrm{C}\mathrm{C}}_i(p_{\mathrm{C}\mathrm{R}})}$ is specific to and generally different for each neuron, reflecting its role in the perceptual decision-making process. A CC stimulus dependence may arise as a result of stimulus-dependent decision feedback (Haefner et al., 2016; Bondy et al., 2018; Lange and Haefner, 2017), or other sources of stimulus-dependent cross-neuronal correlations (Ponce-Alvarez et al., 2013; Orbán et al., 2016) such as shared gain fluctuations (Goris et al., 2014). In fact, we will show below that gain-induced stimulus-dependent cross-neuronal correlations account for observed features in our empirical data. Note that we do not distinguish between CC stimulus dependencies and a dependence of the CC on $p_{\text{CR}}$. We do not make this distinction here because most generally a change in the stimulus level results in a change of $p_{\text{CR}}$, and the two cannot be disentangled. However, the $p_{\text{CR}}$ more generally depends on other factors such as the reward value, attention level, or arousal state, and in Equation 7 the separate dependencies on the stimulus and $p_{\text{CR}}$ can be explicitly indicated as ${\mathrm{CC}}_i(p_{\text{CR}},s)$ when the experimental paradigm allows to separate these two influences.

For simplicity, we presented above only the general relationship between the CP and CC in Equation 7 derived as a linear approximation for weak activity-choice covariations, as this is the regime relevant for single sensory neurons. See Methods for the exact analytical solution from the threshold model (Equation 16) and Appendix 1 for its derivation. Despite the assumption of weak activity-choice covariations, this approximation is very close over the empirically relevant range of CC’s (Figure 2A–B). Below we will focus on a concrete type of CC stimulus dependence, namely originated by gain fluctuations, but it is clear from Equation 7 that any CC stimulus dependence will modify the $\mathrm{CP}(p_{\text{CR}})$ shape induced purely by the threshold effect. A summary of the overall relation between the CP, $\mathrm{CTA}$, and $\mathrm{CC}$ is provided in Figure 2D.

Figure 2

Download asset Open asset

The model provides a concrete prediction of a stereotyped dependence of CP on $p_{\text{CR}}$ through $h(p_{\text{CR}})$ when the choice-related signals are mediated by an intermediate decision variable $d$, which is testable using data. First, under the assumption that CC is constant and therefore $h(p_{\text{CR}})$ is the only source of CP dependence on $p_{\text{CR}}$, for a positive CC ($\mathrm{CP}>0.5$), the $\mathrm{CP}(p_{\text{CR}})$ should have a minimum at $p_{\text{CR}}=0.5$ and increase symmetrically as $p_{\text{CR}}$ deviates from 0.5 as the result of a change in the stimulus in either direction (Figure 2A). When the CC is negative ($\mathrm{CP}<0.5$), then $\mathrm{CP}(p_{\text{CR}})$ should have a maximum at $p_{\text{CR}}=0.5$ and analogously decrease symmetrically as $p_{\text{CR}}$ deviates from 0.5. Second, since the influence of $h(p_{\text{CR}})$ is multiplicative, it creates higher absolute differences in the CP across different stimulus levels for cells with a stronger CP (either larger or smaller than 0.5). Third, the dependence on $h(p_{\text{CR}})$ is weak for a wide range of $p_{\text{CR}}$ values (Figure 2A), making it empirically detectable only when including highly informative stimuli in the analysis to obtain $p_{\text{CR}}$ values very different from 0.5. However, for those $p_{\text{CR}}$ values, CP estimates are less reliable, because only for few trials the choice is expected to be inconsistent with the sensory information, meaning that one of the two distributions $p(r_i|D=1)$ or $p(r_i|D=-1)$ is poorly sampled. This means that to detect the $h(p_{\text{CR}})$ modulation for single cells, many trials would be needed for each value of $p_{\text{CR}}$ to obtain good estimates. Because $h(p_{\text{CR}})$ is common to all cells, averaging $\mathrm{CP}(p_{\text{CR}})$ profiles across cells can also improve the estimation. This averaging may also help to isolate the $h(p_{\text{CR}})$ modulation, assuming that cell-specific ${\mathrm{CP}}_i$ stimulus dependencies introduced through choice correlations ${\mathrm{CC}}_i$ are heterogeneous across cells and average out. We refer to Appendix 1 for a detailed analysis of the statistical power for the detection of $h(p_{\text{CR}})$ as a function of the number of trials and cells used to estimate an average $\mathrm{CP}(p_{\text{CR}})$ profile. We will present below (Section ‘Stimulus dependence of choice-related signals in the responses of MT cell’) evidence for the $h(p_{\text{CR}})$ modulation from a re-analysis of the data in Britten et al., 1996.

We will now focus on a concrete source of stimulus-dependent correlations that leads to a non-constant $\mathrm{CC}(p_{\text{CR}})$, namely the effect of gain fluctuations into the stimulus-response relationship (Goris et al., 2014; Ecker et al., 2014; Kayser et al., 2015; Schölvinck et al., 2015). Goris et al., 2014 showed that 75% of the variability in the responses in monkeys MT cells when presented with drifting gratings could be explained by gain fluctuations. We derive the CP dependencies on $p_{\text{CR}}$ in a feedforward model of decision-making (Shadlen et al., 1996; Haefner et al., 2013) that also models the effect of gain fluctuations in the responses. The feedfoward model considers a population of sensory responses, $\overrightarrow{r}=(r_1,\mathrm{\dots },r_n)$, with tuning functions $\overrightarrow{f}(s)=(f_1(s),\mathrm{\dots },f_n(s))$, responses $r_i=f_i(s)+{\xi }_i$, and a covariance structure $\mathrm{\Sigma }$ of the neuron’s intrinsic variability ${\xi }_i$. The responses are read out into the decision variable with a linear decoder 

where $\overrightarrow{w}$ are the read-out weights. The categorical choice $D$ is made by comparing $d$ to a threshold θ. With this model, the general expression of Equation 7 reduces to

where ${(\mathrm{\Sigma }(s)\overrightarrow{w})}_i=\mathrm{cov}(r_i,d)$ and $\mathrm{var}d={\overrightarrow{w}}^{\top }\mathrm{\Sigma }(s)\overrightarrow{w}$. This expression corresponds to the one derived by Haefner et al., 2013, except for $h(p_{\text{CR}})$ and for the fact that we now explicitly indicate the dependence of the correlation structure $\mathrm{\Sigma }(s)$ on the stimulus. The expression relates the CP magnitude to single-unit properties such as the neurometric sensitivity, as well as to population properties, such as the decoder pooling size and the magnitude of the cross-neuronal correlations, which determine CC (Shadlen et al., 1996; Haefner et al., 2013). In particular, if the decoding weights are optimally tuned to the structure of the covariability $\mathrm{\Sigma }(p_{\text{CR}}=0.5)$ at the decision boundary, this results in a proportionality between ${\mathrm{CP}}_i(p_{\text{CR}}=0.5)$ and the neurometric sensitivity of the cells: ${\mathrm{CP}}_i(p_{\text{CR}}=0.5)\propto f_i^{\prime }/{\sigma }_{r_i}$ (Haefner et al., 2013), as has been experimentally observed (Britten et al., 1996; Parker and Newsome, 1998). While this feedfoward model is generic, we concretely study CC stimulus dependencies induced by the effect of global gain response fluctuations in cross-neuronal correlations. Following Goris et al., 2014 we modeled the responses of cell $i$ in trial $k$ as $f_{ik}(s)=g_k{f}_i(s)$, where g_k is a gain modulation factor shared by the population. We assume that the readout weights $\overrightarrow{w}$ are stimulus-independent. As a consequence, the covariance of population responses $\mathrm{\Sigma }$ has a component due to the gain fluctuations:

 where ${\sigma }_G^2$ is the variance of the gain $g$ and $\overline{\mathrm{\Sigma }}$ is the covariance not associated with the gain, which for simplicity we assume to be stimulus independent. The component of the cross-neuronal covariance matrix $\mathrm{\Sigma }$ induced by gain fluctuations is proportional to the tuning curves ($\propto \overrightarrow{f}(s){\overrightarrow{f}}^T(s)$). A deviation $\mathrm{\Delta }s\equiv s-s_0$ of the stimulus from the uninformative stimulus s₀ produces a change $\mathrm{\Delta }\overrightarrow{f}={\overrightarrow{f}}^{\prime }(s_0)\mathrm{\Delta }s$ in the population firing rates, which affects the variability of the responses, the variability of the decoder, and their covariance, which all vary with $\mathrm{\Delta }s$. Because the variance of the decoder $\mathrm{var}d={\overrightarrow{w}}^{\top }\mathrm{\Sigma }(s)\overrightarrow{w}$ and the covariance $\mathrm{cov}(r_i,d)={(\mathrm{\Sigma }(s)\overrightarrow{w})}_i$ both depend on the concrete form of the read-out weights, the effect of gain-induced stimulus dependencies on the CP is specific for each decoder. Under the assumption of an optimal linear decoder at the decision boundary s₀ (${\displaystyle \overrightarrow{w}\propto {\mathrm{\Sigma }}^{-1}{\overrightarrow{f}}^{\mathrm{\prime }}(s_0)}$), we obtain an approximation of the CC dependence on the stimulus deviation $\mathrm{\Delta }s$ from s₀ (see Methods for details):

where the slope is determined by the coefficient ${\displaystyle {\beta }_{p_{\mathrm{C}\mathrm{R}}}={\sigma }_G{\lambda }_i\left[1-{\mathrm{C}\mathrm{C}}_i^2(p_{\mathrm{C}\mathrm{R}}=0.5)\right]}$, with ${\lambda }_i$ being the fraction of the variance of cell $i$ caused by the gain fluctuations (Methods). The choice rate $p_{\text{CR}}$ is determined by the stimulus $\mathrm{\Delta }s$ as characterized by the psychometric function. For this form of the slope coefficient ${\beta }_{p_{\text{CR}}}$ obtained with an optimal decoder all the factors contributing to it are positive (Figure 2C). In Appendix 4 we further analytically describe how gain fluctuations introduce CP stimulus dependencies not only for an optimal decoder, but also for any unbiased decoders. Conversely to the factor $h(p_{\text{CR}})$, the pattern of $\mathrm{CP}(p_{\text{CR}})$ profiles produced by the gain fluctuations is cell-specific, with a stronger asymmetric component for cells with higher ${\lambda }_i$ (Figure 2C). Furthermore, while the sign of the multiplicative modulation $h(p_{\text{CR}})$ changes when CC>0 or CC<0, the gain-induced contribution in Equation 11 is additive. As seen in Figure 2C, for cells with a weak activity-choice covariation for uninformative stimuli ($\mathrm{CP}(p_{\text{CR}}=0.5)$ close to 0.5), this implies that the CP of a neuron can actually change from below 0.5 to above 0.5 across the stimulus range presented in the experiment.

In the light of our findings above, we re-analyzed the classic Britten et al., 1996 data containing responses of neurons in area MT in a coarse motion direction discrimination task (see Methods for a description of the data set). Our objective is to identify any patterns of CP dependence on the choice rate/stimulus level. First, we describe our results testing for the threshold-induced CP stimulus dependence, $h(p_{\text{CR}})$, and then more generally we characterize the $\mathrm{CP}(p_{\text{CR}})$ patterns found in the data using clustering analysis. Finally, as an alternative to CP analysis, we show how to extend Generalized Linear Models (GLMs) of neural activity to include stimulus-choice interaction terms that incorporate the stimulus dependencies of activity-choice covariations derived with our theoretical approach and found above in the MT data.

Testing the presence of a threshold-induced CP stimulus dependence in experimental data

We start describing how to analyze within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles to test the existence of the threshold-induced modulation. The theoretically derived properties of $h(p_{\text{CR}})$ suggest several empirical signatures that will be reflected in the within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles. First, because $h(p_{\text{CR}})$ introduces a multiplicative modulation of the choice correlation, for informative stimuli it leads to an increase of the CP for cells with positive choice correlation ($\mathrm{CP}>0.5$) and to a decrease for cells with negative choice correlation ($\mathrm{CP}<0.5$). Second, because $h(p_{\text{CR}})$ is multiplicative, the absolute magnitude of the modulation will be higher for cells with stronger choice correlation, that is CPs most different from 0.5. Third, the effect of $h(p_{\text{CR}})$ is strongest when one choice dominates and hence most noticeable for highly informative stimuli.

These properties of $h(p_{\text{CR}})$ indicate that, to detect this modulation, it is necessary to examine within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles isolated from across-cells heterogeneity in the magnitude of the CP values. Ideally, we would like to calculate a $\mathrm{CP}(p_{\text{CR}})$ profile for each cell and analyze the shape of these single-cell profiles. However, given the available number of trials, estimates of ${\displaystyle \mathrm{C}\mathrm{P}(p_{\mathrm{C}\mathrm{R}})}$ profiles for single cells are expected to be noisy. The estimation error of the CP is higher when $p_{\text{CR}}$ is close to 0 or 1, the same $p_{\text{CR}}$ values for which the ${\displaystyle h(p_{\mathrm{C}\mathrm{R}})}$ modulation would be most noticeable. The standard error of $\widehat{\mathrm{CP}}$ can be approximated as ${\displaystyle \mathrm{S}\mathrm{E}\mathrm{M}(\widehat{\mathrm{C}\mathrm{P}})\approx 1/\sqrt{12K{p}_{\mathrm{C}\mathrm{R}}(1-p_{\mathrm{C}\mathrm{R}})}}$ (Bamber, 1975; Hanley and McNeil, 1982, see Methods), where $K$ is the number of trials. In the Britten et al. data set the number of trials varies for different stimulus levels, and most frequently $K=30$ for highly informative stimuli. In that case, for $p_{\text{CR}}=0.9$, only three trials for choice $D=-1$ are expected, and $\mathrm{SEM}(\widehat{\mathrm{CP}})\approx 0.18$. As can be seen from Figure 2A, this error surpasses the order of magnitude of the CP modulations expected from $h(p_{\text{CR}})$. This means that we need to combine CP estimates of adjacent $p_{\text{CR}}$ values, and/or combine estimated $\mathrm{CP}(p_{\text{CR}})$ profiles across neurons, to reduce the standard error (See Appendix 1 for a detailed analysis of the statistical power for the detection of $h(p_{\text{CR}})$).

When averaging CPs across neurons, two considerations are important. First, cells that for $p_{\text{CR}}=0.5$ have a CP higher or lower than 0.5 should be separated, given that the sign of the CC leads to an inversion of the profile resulting from $h(p_{\text{CR}})$ (Equation 7). If not separated, the $h(p_{\text{CR}})$-dependence would average out, or the average $\mathrm{CP}(p_{\text{CR}})$ profile would reflect the proportion of cells with CPs higher or lower than 0.5 in the data set. Second, the average should correspond to an average -across cells- of within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles, and hence it should only include cells for which a full $\mathrm{CP}(p_{\text{CR}})$ profile can be calculated. This is important because for each cell $i$ the $h(p_{\text{CR}})$ modulation is relative to the value of ${\mathrm{CP}}_i(p_{\text{CR}}=0.5)$. If a different subset of cells was included in the average of the CP at each $p_{\text{CR}}$ value, the resulting shape across $p_{\text{CR}}$ values of the averaged CPs would not be an average of within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles. Conversely, in that case, the resulting shape would reflect the heterogeneity in the magnitude of the CP values across the subsets of cells averaged at each $p_{\text{CR}}$ value. In the single-cell recordings from Britten et al., the range of stimulus levels used varies across neurons, and for a substantial part of the cells a full $\mathrm{CP}(p_{\text{CR}})$ profile cannot be constructed. Following the second consideration, those cells were excluded from the analysis to avoid that they only contributed to the average at certain $p_{\text{CR}}$ values.

We derived the following refined procedure to analyze $\mathrm{CP}(p_{\text{CR}})$ profiles. As a first step, we constructed a $\mathrm{CP}(p_{\text{CR}})$ profile for each cell. First, for each cell and each stimulus coherence level we calculated a CP estimate if at least four trials were available for each decision. For the experimental data set, CPs are always estimated from its definition (Equation 1), and we will only use the theoretical expression of $h(p_{\text{CR}})$ to fit the modulation of the experimentally estimated $\mathrm{CP}(p_{\text{CR}})$ profiles. Second, as a first way to improve the CP estimates, we binned $p_{\text{CR}}$ values into five bins and assigned stimulus coherence levels to the bins according to the psychometric function that maps stimulus levels to $p_{\text{CR}}$, with the central bin containing the trials from the zero-signal stimulus. A single CP value per bin for each cell was then obtained as a weighted average of the CPs from stimulus levels assigned to each bin. The weights were calculated as inversely proportional to the standard error of the estimates, giving more weight to the most reliable CPs (see Methods). The results that we present hereafter are all robust to the selection of the minimum number of trials and the binning intervals. Unless otherwise stated, in all following analyses we included all the cells ($N=107$) for which we had data to compute CPs in all five bins, thus allowing us to estimate a full within-cell $\mathrm{CP}(p_{\text{CR}})$ profile. As a second step, we averaged the within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles across cells, taking into account the two considerations above. As before, averages were weighted by inverse estimation errors.

Figure 3A shows the averaged $\mathrm{CP}(p_{\text{CR}})$ profiles. To assess the statistical significance of the CP dependence on $p_{\text{CR}}$, we developed a surrogates method to test whether a pattern consistent with the predicted CP-increase for informative stimuli could appear under the null hypothesis that the CP has a constant value independent of $p_{\text{CR}}$ (see Methods). For the cells with average CP higher than 0.5, we found that the modulation of the CP was significant ($p=0.0006$), with higher CPs obtained for $p_{\text{CR}}$ close to 0 or one in agreement with the model. For cells with average CP lower than 0.5, the modulation was not significant ($p=0.26$). While the actual absence of a modulation would imply that the choice-related signals in these neurons are not mediated by a continuous intermediate decision-variable but may be, for example, due to categorical feedback, we point out the lower power of this statistical test due to fewer neurons being in the $\mathrm{CP}<0.5$ group and the expected effect size being lower, too. First, there were 74 cells with CP higher than 0.5 but only 33 with CP lower than 0.5, meaning that the estimation error is larger for the average $\mathrm{CP}(p_{\text{CR}})$ profile of the cells with $\mathrm{CP}<0.5$. Second, as the modulation predicted by $h(p_{\text{CR}})$ is multiplicative, its impact is expected to be smaller when the magnitude of $\mathrm{CP}-0.5$ is smaller. Figure 3A shows that CP values are on average closer to 0.5 for the cells with $\mathrm{CP}<0.5$, in agreement with Figure 5 of Britten et al., 1996. This means that fewer cells classified in the group with $\mathrm{CP}<0.5$ have choice-related responses. Therefore, the fact that we cannot validate the prediction of an inverted symmetric $h(p_{\text{CR}})$ modulation for the cells with $\mathrm{CP}<0.5$ with respect to the cells with $\mathrm{CP}>0.5$ is not strong evidence against the existence of a threshold-induced CP stimulus dependence. We further confirmed the robustness of the results in a wider set of cells. For this purpose, we repeated the analysis forming subsets separately including cells with a computable CP for the three bins with $p_{\text{CR}}$ lower or equal 0.5, and the three with $p_{\text{CR}}$ higher or equal than 0.5. Also in this case the observed $\mathrm{CP}(p_{\text{CR}})$ pattern was significant ($p=0.0013$) for cells with average CP higher than 0.5 (Figure 3B, $N=171$), and non-significant for cells with CP lower than 0.5 (p=0.20).

Figure 3

Download asset Open asset

Interestingly, the identified significant $\mathrm{CP}(p_{\text{CR}})$ dependence for the cells with $\mathrm{CP}>0.5$ goes beyond the symmetric threshold-induced shape predicted by $h(p_{\text{CR}})$, both in magnitude and shape (Figure 2A), since the increase is bigger for $p_{\text{CR}}$ values close to 1 than to 0. This implies that the choice correlation for each neuron, ${\mathrm{CC}}_i(p_{\text{CR}})$, must systematically change with $p_{\text{CR}}$ as well, contributing to the overall CP stimulus dependency observed. In particular, the observed average $\mathrm{CP}(p_{\text{CR}})$ profile indicates that the CP increase appears to be higher for $p_{\text{CR}}>0.5$. The finding of this asymmetry is consistent with results reported in Britten et al., 1996, who found a significant but modest effect of coherence direction on the CP (see their Figure 3). By experimental design, the direction of the dots corresponding to choice $D=1$ was tuned for each cell separately to coincide with their most responsive direction. This means that this asymmetry indicates that CPs tend to increase more when the stimulus provides evidence for the direction eliciting a higher response. However, Britten et al., 1996 found no significant relation between the global magnitude of the firing rate and the CP (see their Figure 3), and we confirmed this lack of relation specifically for the subset of $N=107$ cells (no significant correlation coefficient between average rate and average CP values, $p=0.33$). This eliminates the possibility that higher CPs for high $p_{\text{CR}}>0.5$ values are due only to higher responses, and suggests a richer underlying structure of $\mathrm{CP}(p_{\text{CR}})$ patterns, which we will investigate next using cluster analysis to identify the predominant patterns shared by the within-cell $\mathrm{CP}(p_{\text{CR}})$ profiles.

Characterizing the experimental patterns of CP stimulus dependencies with cluster analysis

We carried out unsupervised $k$-means clustering (Bishop, 2006) to examine the patterns of $\mathrm{CP}(p_{\text{CR}})$ without a priori assumptions about a modulation $h(p_{\text{CR}})$ associated with the threshold effect. Clustering was performed on $\mathrm{CP}(p_{\text{CR}})-0.5$, with each cell represented as a vector in a five-dimensional space, where five is the number of $p_{\text{CR}}$ bins used to summarize the data as described above. To consider both the shape and sign of the modulation, distances between neurons were calculated with the cosine distance between their $\mathrm{CP}(p_{\text{CR}})$ profiles (one minus the cosine of the angle between the two vectors). Clustering was performed for a range of specified numbers of clusters. Specifying the existence of two clusters, we naturally recovered the distinction between cells with CP higher or lower than 0.5 (Figure 4A). The statistical significance of any ${\displaystyle p_{\mathrm{C}\mathrm{R}}}$-modulation was again assessed constructing surrogate $\mathrm{CP}(p_{\text{CR}})$ profiles and repeating the clustering analysis on those surrogates. As before, a significant dependence of the CP on $p_{\text{CR}}$ was found only for the cluster associated with CP higher than 0.5 ($p=0.0007$ for $\mathrm{CP}>0.5$ and $p=0.21$ for $\mathrm{CP}<0.5$).

Figure 4

Download asset Open asset

As mentioned above, the divergence from $h(p_{\text{CR}})$ of the average ${\displaystyle \mathrm{C}\mathrm{P}(p_{\mathrm{C}\mathrm{R}})}$ profile for cells with $\mathrm{CP}>0.5$ suggests that cell-specific modulations are introduced through ${\mathrm{CC}}_i(p_{\text{CR}})$. While the variability of individual ${\mathrm{CP}}_i(p_{\text{CR}})$ profiles (Figure 3C) is expected to reflect substantially the high estimation errors of $\widehat{\mathrm{CP}}$ for the single cells (Figure 3D), the presence of subclusters can identify $\mathrm{CP}(p_{\text{CR}})$ patterns common across cells.

We proceed to examine subclusters within the $\mathrm{CP}>0.5$ cluster with a significant ${\displaystyle \mathrm{C}\mathrm{P}(p_{\mathrm{C}\mathrm{R}})}$ profile, excluding from our analysis cells within the $\mathrm{CP}<0.5$ cluster (analogous results were found when increasing the number of clusters in a nonhierarchical way, without a priori excluding these cells, see Appendix 3—figure 1A). Average $\mathrm{CP}(p_{\text{CR}})$ profiles obtained when inferring two subclusters of cells with $\mathrm{CP}>0.5$ are shown in Figure 4B. For both subclusters the $\mathrm{CP}(p_{\text{CR}})$ dependence is significant ($p=0.0008$ for cluster two and $p=0.0026$ for cluster 3, respectively, in Figure 4B). The larger cluster has a more symmetric shape of dependence on $p_{\text{CR}}$, with an increase of CP in both directions when the stimulus is informative, consistent with the prediction of a threshold-induced CP stimulus dependence $h(p_{\text{CR}})$. For the smaller cluster the dependence is asymmetric, with a CP increase when the stimulus direction is consistent with the preferred direction of the cells and a decrease in the opposite direction. We verified that no significant difference exists between the firing rates of the cells in the two subclusters (Wilcoxon rank-sum test, $p=0.23$). The monotonic shape of the second subcluster mirrors the dependency produced by response gain fluctuations as predicted by the gain model described above. This suggests that the neurons in this subcluster differ from the neurons in the other subcluster by a substantially larger gain-induced variability, a testable prediction for future experiments and further discussed below.

Introducing a second cluster allows for representing each neuron’s $\mathrm{CP}(p_{\text{CR}})$-dependency in the two-dimensional space (Figure 4C) spanned by the mean profiles for each of the three clusters. The horizontal axis corresponds to the separation between the two initial clusters, and is closely aligned to the departure of the average CP from 0.5. The vertical axis is defined by the vectors corresponding to the centers of the two subclusters and hence is determined separately for the cells with average CP higher and lower than 0.5 (see Methods for details, and Appendix 3—figure 1A). The vertical axis is associated with the degree to which the $\mathrm{CP}(p_{\text{CR}})$ dependence is symmetric or asymmetric with respect to $p_{\text{CR}}=0.5$. Cells for which the CP increases consistently with its preferred direction of motion coherence lie on the upper half-plane. To further support this interpretation of the axis, we repeated the clustering procedure replacing the nonparametric $k$-means procedure with a parametric procedure that defines the subclusters with a symmetric and an asymmetric template, respectively. The data is distributed approximately equally in both spaces (Figure 4C–D).

Similar results were also obtained when increasing the number of clusters non-hierarchically. Introducing a third cluster for all cells leaves almost unaltered the cluster of cells with CP lower than 0.5 (Appendix 3—figure 1B). The cluster of cells with CP higher than 0.5 splits into two subclusters analogous to the ones found from cells with CP higher than 0.5 alone. The distinction between cells with more symmetric and asymmetric $\mathrm{CP}(p_{\text{CR}})$ dependencies is robust to the selection of a larger number of clusters, that is, clusters with this type of dependencies remain large when allowing for the discrimination of more patterns (Appendix 3—figure 1C). However, we do not mean to claim that the variety of $\mathrm{CP}(p_{\text{CR}})$ profiles across cells can be reduced to three separable clusters. As reflected in the distributions in Figure 4C–D, the clusters are not neatly separable. Indeed, a richer variety of profiles would be expected if the properties of $\mathrm{CP}(p_{\text{CR}})$ profiles across cells were associated with their tuning properties and the structure of feedback projections, as we further argue in the Discussion. The predominance of a symmetric and asymmetric pattern would only reflect which are the predominant $\mathrm{CP}(p_{\text{CR}})$ shapes shared across cells.

This clustering analysis confirms the presence of shared patterns of CP stimulus-dependence across cells, whose shape is compatible with the analytical predictions from the threshold- and gain-related dependencies. The symmetric component of CP stimulus dependence is congruent with $h(p_{\text{CR}})$ (Equation 6), albeit with a larger magnitude than predicted (Figures 2A and 3A, and additional analysis of the statistical power in Appendix 1). This stronger modulation suggests an additional symmetric contribution of the choice correlation $\mathrm{CC}(p_{\text{CR}})$ and/or a dynamic feedback reinforcing the stronger modulation for highly informative stimuli. However, while the cluster analysis separates the predominant $\mathrm{CP}(p_{\text{CR}})$ patterns, the Britten et al. data lacks the statistical power to further distinguish between $h(p_{\text{CR}})$ and symmetric $\mathrm{CC}(p_{\text{CR}})$ contributions with a similar shape.

The implications of a stimulus-dependent relationship between the behavioral choice and sensory neural responses are not restricted to measuring them as CPs, for which activity-choice covariations are quantified without incorporating other explanatory factors of neural responses. To further substantiate the existence of this stimulus-dependent relationship in MT data, and to understand how our model predictions could help to refine other analytical approaches, we examined how representing that relationship can improve statistical models of neural responses. In particular, we study how the stimulus-dependent choice-related signals that we discovered may inform the refinement of Generalized Linear Models (GLMs) of neural responses (Truccolo et al., 2005; Pillow et al., 2008). In the last few years, GLMs have been used for modelling choice dependencies together with the dependence on other explanatory variables, such as the external stimulus, response memory, or interactions across neurons (Park et al., 2014; Runyan et al., 2017). Typically, in a GLM of firing rates each explanatory variable contributes with a multiplicative factor that modulates the mean of a Poisson process. In their classical implementation, the choice modulates the firing rate as a binary gain factor, with a different gain for each of the two choices (Park et al., 2014; Runyan et al., 2017; Pinto et al., 2019). The multiplicative nature of this factor already introduces some covariation between the impact of the choice on the rate and the one of the other explanatory variables. However, using a single regression coefficient to model the effect of the choice on the neural responses may be insufficient if choice-related signals are stimulus dependent, as suggested by our theoretical and experimental analysis.

We developed a GLM (see Methods) that can model stimulus-dependencies of choice signals (or, in other words, stimulus-choice interactions) by including multiple choice-related predictors that allow for a different strength of dependence of the firing rate on the choice for different subsets of stimulus levels (via the choice rate, $p_{\text{CR}}$). We fitted this model, which we call the stimulus-dependent-choice GLM, to MT data and we compared its cross-validated performance against two traditional GLMs. In the first type, called the stimulus-only GLM, the rate in each trial is predicted only based on the external stimulus level. In a second type, that we called stimulus-independent-choice GLM and that corresponds to the traditional way to include choice signals in a GLM (Park et al., 2014; Runyan et al., 2017; Scott et al., 2017; Pinto et al., 2019; Minderer et al., 2019), additionally the effect of choice is included, but using only a single, stimulus-independent choice predictor.

To compare the models, we separated the trials recorded from each MT cell (Britten et al., 1996) into training and testing sets, and calculated the average cross-validated likelihood for each type of model on the held-out testing set. To quantify the increase in predictability when adding the choice as a predictor we defined the relative increase in likelihood (RIL) as the relative increase of further adding the choice as a predictor relative to the increase of previously adding the stimulus as a predictor. RIL measures the relative influence of the choice and the sensory input in the neural responses. Figure 5A compares the cross-validated RIL values obtained on MT neural data when fitting either the stimulus-independent-choice or the stimulus-dependent-choice GLMs. We found that RIL values were mostly higher when allowing for multiple choice parameters, both in terms of average RIL values (Figure 5C) and in terms of the proportion of cells in each cluster for which the RIL was higher than a certain threshold, here selected to be at 10% (Figure 5B).

Figure 5

Download asset Open asset

GLMs that include stimulus-choice interaction terms can be used not only to better describe the firing rate of neural responses, but also to individuate more precisely the neurons or areas by their choice signals. To illustrate this point, we show how adding the interaction terms may change the relative comparison of cells by their RIL values. Consider the three neurons with highest RIL for the stimulus-dependent-choice GLM (Figure 5A, and with corresponding $\mathrm{CP}(p_{\text{CR}})$ profiles shown in Figure 5D). The ranking of cells 1 and 2 by RIL flips with respect to the stimulus-independent-choice GLM because of the higher $\mathrm{CP}(p_{\text{CR}})$ modulation of cell 2. Similarly, while the RIL with multiple choice parameters for cells 1 and 3 are close, the RIL of cell 3 is substantially lower with a single choice parameter, indicating that its pattern of stimulus dependence is less well captured by a single parameter. The degree to which a model with interaction terms improves the predictability will depend on the shape of the $\mathrm{CP}(p_{\text{CR}})$ patterns, which themselves are expected to vary across areas or across cells with different tuning properties. For example, we see in Figure 5C that for the cluster with an asymmetric $\mathrm{CP}(p_{\text{CR}})$ profile (cluster 3), the average RIL with only one choice parameter suggests that this type of cells are not choice driven. The reason is that for the cells in this cluster the sign of the choice influence on the rate can be stimulus dependent, which is impossible to model by a single choice parameter. Furthermore, the profile of the GLM choice parameters across stimulus levels provides a characterization of stimulus-dependent choice-related signals analogous to the $\mathrm{CP}(p_{\text{CR}})$ profile, in this case within the GLM framework, hence allowing efficient inference including principled regularization and the ability to account for a range of factors beyond choices and stimuli. Overall, we expect that accounting for stimulus-choice interactions in GLMs will allow for a more accurate assessment of the relative importance of stimulus and choice on neural responses.

Previous work has demonstrated that solving analytically models of perceptual decision-making can lead to important new insights on the interpretation of the relationship between neural activity and choice in terms of decision-making computations (Bogacz et al., 2006; Gold and Shadlen, 2007; Haefner et al., 2013). In particular, previous analytical work on CPs has shown how experimentally measured CPs relate to the read-out weights by which sensory neurons contribute to the internal stimulus decoder in a feedforward model, assuming both choices are equally likely (Haefner et al., 2013; Pitkow et al., 2015). Here, we provided a general analytical solution of CPs in a more general model, with informative stimuli resulting in an unbalanced choice rate, and valid both for feedforward and feedback choice signals. We derived the analytical dependency of CP on the probability of one of the choices ($p_{\text{CR}}\equiv p(\mathrm{choice}=1)$), which mediates the dependence of the CP on the stimulus strength. Our model is therefore directly applicable to both discrimination and detection tasks, for any stimulus strength that elicits both choices. As we demonstrated, these advances in the analytical solution of the decision-threshold model allowed for detecting and interpreting stimulus dependencies of choice-related signals in neural activity.

We first derive our analytical CP expression valid in the presence of informative stimuli, decision-related feedback, and top-down sources of activity-choice covariation, such as prior bias, trial-to-trial memory, or internal state fluctuations. We follow Haefner et al., 2013 and assume a threshold model of decision making, in which the choice $D$ is triggered by comparing a decision variable $d$ with a threshold θ, so that if $d>\theta$ choice $D=1$ is made, and $D=-1$ otherwise. The identification of the binary choices as $D=\pm 1$ is arbitrary and an analogous expression would hold with another mapping of the categorical variable. The choice probability (Britten et al., 1996) of cell $i$ is defined as

and measures the separation between the two choice-specific response distributions $p(r_i|D=-1)$ and $p(r_i|D=1)$. It quantifies the probability of responses to choice $D=1$ to be higher than responses to $D=-1$. If there is no dependence between the choice and the responses this probability is $\mathrm{CP}=0.5$. To obtain an exact solution of the CP, we assume that the distribution $p(r_i,d)$ of the responses r_i of cell $i$ and the decision variable $d$ can be well approximated by a bivariate Gaussian. Under this assumption, following Haefner et al., 2013 (see their Supplementary Material) the probability of the responses for choice $D=1$ follows the distribution 

where a more parsimonious expression is obtained using the z-score $z_i=(r_i-⟨r_i⟩)/{\sigma }_{r_i}$. This distribution is a skew-normal (Azzalini, 1985), where $\varphi (\cdot ;0,1)$ is the standard normal distribution with zero mean and unit variance, and $\mathrm{\Phi }(\cdot ;0,1)$ is its cumulative function. Furthermore, ${\sigma }_{r_i,d}$ is the covariance of r_i and $d$, ${\sigma }_{d|r_i}$ is the conditional standard deviation of $d$ given r_i, and the probability of $D=1$ is

which determines the rate of each choice over trials. The choice $D=-1$ could equally be taken as the choice of reference, resulting in an analogous formulation. Intuitively, $p_{\text{CR}}$ increases when the mean of the decision variable $⟨d⟩$ is higher than the threshold θ, and decreases when its standard deviation ${\sigma }_d$ increases. The form of the distribution of Equation 13 can be synthesized in terms of $p_{\text{CR}}$ and the correlation coefficient ${\rho }_{r_{i}d}$, which was named by Pitkow et al., 2015 choice correlation (${\mathrm{CC}}_i$). In particular, defining $\alpha \equiv {\rho }_{r_{i}d}/\sqrt{1-{\rho }_{r_{i}d}^2}$ and $c\equiv {\mathrm{\Phi }}^{-1}(p_{\text{CR}})/\sqrt{1-{\rho }_{r_{i}d}^2}$

The CP is completely determined by $p(z_i|D=-1)$ and $p(z_i|D=1)$, and these distributions depend only on $p_{\text{CR}}$ and the correlation coefficient ${\rho }_{r_{i}d}$. Plugging the distribution of Equation 15 into the definition of the CP (Equation 12) an analytical solution is obtained:

where $\mathrm{T}$ is the Owen’s T function (Owen, 1956). In Appendix 1, we provide further details of how this expression is derived. For an uninformative stimulus ($p_{\text{CR}}=0.5$), the function $\mathrm{T}$ reduces to the arctangent and the exact result obtained in Haefner et al., 2013 is recovered. The dependence on ${\rho }_{r_{i}d}$ can be understood because under the Gaussian assumption the linear correlation captures all the dependence between the responses and the decision variable $d$. The dependence on $p_{\text{CR}}$ reflects the influence of the threshold mechanism, which maps the dependence of r_i with $d$ into a dependence with choice $D$ by partitioning the space of $d$ in two regions.

While Equation 16 provides an exact solution of the CP, in the Results section we present and mostly focus on a linear approximation to understand how the stimulus content modulates the choice probability. This approximation is derived (Appendix 1) in the limit of a small ${\rho }_{r_{i}d}$, which leads to CPs close to 0.5 as usually measured in sensory areas (Nienborg et al., 2012). However, as we show in the Results and further justify in the Appendix this approximation is robust for a wide range of ${\rho }_{r_{i}d}$ values. The linear approximation relates the choice probability to the Choice Triggered Average (CTA) (Haefner, 2015; Chicharro et al., 2017), defined as the difference of the mean responses for each choice (Equation 3). Given the binary nature of choice $D$, the $\mathrm{CTA}$ is directly proportional to the covariance of the responses and the choice: ${\mathrm{CTA}}_i=\mathrm{cov}(r_i,D)/[2p(D=1)p(D=-1)]$. [Note: This relation holds for the covariance between any variable $x$ and a binary variable $D$, and independently of the convention adopted for the values of $D$: the factor 2 has to be replaced by $a-b$ in general for $D=a,b$ instead of $D=1,-1$.] This relation between ${\mathrm{CTA}}_i$ and $\mathrm{cov}(r_i,D)$, given the factorization $\mathrm{corr}(r_i,D)={\mathrm{CC}}_i\mathrm{corr}(d,D)$ resulting from the mediating decision variable $d$ in the threshold model, allows expressing the ${\mathrm{CTA}}_i$ as in Equation 5, connecting ${\mathrm{CTA}}_i$ to the choice-triggered average of $d$, ${\mathrm{CTA}}_d$. This connection indicates that in the threshold model ${\mathrm{CTA}}_i$ is expected to be stimulus dependent, since an informative stimulus $s$ shifts the mean of $d$, thus altering the dichotomization of $d$ produced by the threshold θ. The exact form of ${\mathrm{CTA}}_d$ depends on the distribution $p(d)$. However, since $d$ is determined by a whole population of neurons, its distribution is expected to be well approximated by a Gaussian distribution, even if the distribution of neural responses for any single neuron is not Gaussian. With this Gaussian approximation, the normalized ${\mathrm{CTA}}_d$ in Equation 5, namely ${\mathrm{CTA}}_d/\sqrt{\mathrm{var}d}$, is specified in terms of the probability of choosing choice 1, $p_{\text{CR}}\equiv p(D=1)=p(d>\theta )$, by the factor $h(p_{\text{CR}})$ (Equation 6). In more detail, the CTA is

To study stimulus dependencies in the relationship between the responses of sensory neurons and the behavioral choice, we analyzed the data from Britten et al., 1996 publicly available in the Neural Signal Archive (http://www.neuralsignal.org). In particular, we analyzed data from file nsa2004.1, which contains single unit responses of macaque MT cells during a random dot discrimination task. This file contains 213 cells from three monkeys. We also used file nsa2004.2, which contains paired single units recordings from 38 sites from one monkey. In the experimental design, for the single unit recordings the direction tuning curve of each neuron was used to assign a preferred-null axis of stimulus motion, such that opposite directions along the axis yield a maximal difference in responsiveness (Bair et al., 2001). For paired recordings, the direction of stimulus motion was selected based on the direction tuning curve of the two neurons and the criterion used to assign it varied depending on the similarity between the tuning curves. For cells with similar tuning, a compromise between the preferred directions of the two neurons was made. For cells with different tuning, the axis were chosen to match the preference of the most responsive cell. To minimize the influence in our analysis of the direction of motion selection, we only analyzed the most responsive cell from each site. Accordingly, our initial data set consisted in a total of 251 cells. The same qualitative results were obtained when limiting the analysis to data from nsa2004.1 alone. Further criteria regarding the number of trials per each stimulus level were used to select the cells. As discussed below, if not indicated otherwise, we present the results from 107 cells that fulfilled all the criteria required.

Our analysis of choice probabilities stimulus dependencies is based on examining the patterns in the $\mathrm{CP}(p_{\text{CR}})$ profile as a function of the probability $p_{\text{CR}}\equiv p(D=1)$. We here describe how these profiles are constructed, the surrogates-based method used to assess the significance of stimulus dependencies, and the clustering analysis used to identify different stimulus dependence patterns. Matlab functions are available at https://github.com/DanielChicharro/CP_DP (Chicharro, 2021; copy archived at swh:1:rev:5850c573860eb04317e7dc550f96b1f47ca91c6a) to calculate weighted average CPs, to obtain CP profiles, and to generate surrogates consistent with the null hypothesis of a constant CP.

To model the effect on the CP of response gain fluctuations we adopted a classic feedforward encoding/decoding model (Shadlen et al., 1996; Haefner et al., 2013), with a linear decoder $d={\overrightarrow{w}}^{\top }\overrightarrow{r}$ (Equation 8), for which the CP depends on cross-neuronal correlations and the read-out weights $\overrightarrow{w}$ following Equation 9. This expression can be derived from Equation 7 directly calculating the choice correlation from its definition (Equation 2). The expressions $\mathrm{cov}(r_i,d)={(\mathrm{\Sigma }(s)\overrightarrow{w})}_i$ and ${\sigma }_d^2={\overrightarrow{w}}^{\top }\mathrm{\Sigma }(s)\overrightarrow{w}$ are obtained as derived in the Supplementary Material S1 of Haefner et al., 2013. For this model, if the read-out weights are optimized to the form of covariability for the uninformative stimulus s₀ at the decision boundary, the CPs are proportional to the neurometric sensitivity of the cells (Haefner et al., 2013; Pitkow et al., 2015), a relationship for which there is some experimental support (e.g. Britten et al., 1996; Parker et al., 2002, reviewed in Nienborg et al., 2012). In more detail, modeling the responses as $r_i=f_i(s)+{\xi }_i$, with tuning functions $\overrightarrow{f}(s)=(f_1(s),\mathrm{\dots },f_n(s))$ and a covariance structure $\mathrm{\Sigma }$ of the neuron’s intrinsic variability ${\xi }_i$, the optimal read-out weights have the form

where ${\overrightarrow{f}}^{\prime }(s_0)$ and $\mathrm{\Sigma }(s_0)$ are the derivative of the tuning curves and the responses covariance matrix, respectively, for $s=s_0$. With these optimal weights, the covariability of the population responses is unbundled, with ${\mathrm{\Sigma }}^{-1}(s_0)$ canceling the effect of $\mathrm{\Sigma }(s_0)$ in $\mathrm{cov}(r_i,d)={(\mathrm{\Sigma }(s)\overrightarrow{w})}_i$, and for each cell the CC is proportional to its own neurometric sensitivity $f_i^{\prime }/{\sigma }_{r_i}$, namely

While this expression is valid for the uninformative stimulus s₀, we examined how this CP expression is perturbed for other informative stimuli $s$ in the presence of gain fluctuations that make the covariance structure $\mathrm{\Sigma }(s)$ stimulus-dependent, altering the structure for which the read-out weights are optimized. Goris et al., 2014 estimated that in MT gain fluctuations accounted for more than 75% of the variance in the responses to sinusoidal gratings, and we found that in the data set of Britten et al., 1996 gain fluctuations also explain a large fraction of trial-to-trial variability of the neurons (62 ± 25% across neurons). Trial-to-trial excitability fluctuations are modeled as a a gain modulatory factor g_k, such that the tuning function for cell $i$ in trial $k$, is $f_{ik}(s)=g_k{f}_i(s)$. In general, the magnitude of the gain may vary across cells, as well as the degree to which the gain co-fluctuates across cells. We here modeled a global gain fluctuation affecting the response of the whole population. Given the gain variability, the covariance structure can be partitioned as in Equation 10, as the sum of a component $\overline{\mathrm{\Sigma }}$ unrelated to the gain fluctuations – which for simplicity we consider to be stimulus-independent – and the gain-induced covariance ${\sigma }_G^2\overrightarrow{f}(s){\overrightarrow{f}}^T(s)$. In a first order approximation, a change $\mathrm{\Delta }s=s-s_0$ in the stimulus leads to a change in the covariance structure such that

where $\mathrm{\Sigma }(s_0)=\overline{\mathrm{\Sigma }}+{\sigma }_G^2\overrightarrow{f}(s_0){\overrightarrow{f}}^T(s_0)$ is the covariance structure for which the weights are optimized. Combining this covariance structure with the form of the optimal read-out weights (Equation 20), we derive the changes in $\mathrm{cov}(r_i,d)$, ${\sigma }_{r_i}^2$, and ${\sigma }_d^2$ with $\mathrm{\Delta }s$, and given Equation 9 determine the perturbation of the CP, leading to the following CP expressions

where ${\sigma }_d^2(s_0)={\overrightarrow{w}}^{\top }\mathrm{\Sigma }(s_0)\overrightarrow{w}$, ${\sigma }_{r_i}^2(s_0)={\mathrm{\Sigma }}_{ii}(s_0)$, ${\beta }_{{\sigma }_r}\equiv {\sigma }_G^2{f}_i(s_0)f_i^{\prime }(s_0)/{\sigma }_{r_i}^2(s_0)$, and ${\lambda }_i^2\equiv {\sigma }_{r_{i}g}^2(s_0)/{\sigma }_{r_i}^2(s_0)$, as introduced in Equation 11, with s₀ resulting in $p_{\text{CR}}=0.5$ for an unbiased decoder. Equation 23a relates the choice correlation for $p_{\text{CR}}=0.5$ with the choice correlation ${\mathrm{CC}}_{i0}(p_{\text{CR}}=0.5)$ that cell $i$ would have if there were no gain fluctuations (${\sigma }_G^2=0$). The coefficient ${\beta }_{\mathrm{CC}}\equiv \sqrt{1-{\lambda }_i^2}$ indicates a decrease in the CC in the presence of gain fluctuations, because of the increase in the response variability produced by the fluctuations, namely ${\sigma }_{r_{i}g}^2={\sigma }_G^2{f}_i^2$. Equation 23b describes the $\mathrm{CC}(p_{\text{CR}})$ profile induced by the gain fluctuations. The second summand corresponds to the increase in the choice correlation due to a new component of $\mathrm{cov}(r_i,d)$ proportional to $\mathrm{\Delta }s$, given that the whole population response determining $d$ is jointly modulated by the gain. In the first summand, the factor $\left[1-{\beta }_{{\sigma }_r}\mathrm{\Delta }s(p_{\text{CR}})\right]$ indicates an attenuation of ${\mathrm{CC}}_i(p_{\text{CR}}=0.5)$ analogous to $\sqrt{1-{\lambda }_i^2}$ in Equation 23a, associated with an increase of the variance in the responses r_i due to $\mathrm{\Delta }s$. Rearranging the terms in Equation 23b, and taking into account the form of the CC for $p_{\text{CR}}=0.5$ (Equation 21), the expression in Equation 11 is obtained, which indicates that the overall effect of the gain fluctuations is an increase of the choice correlation for the stimuli to which the cell is more responsive. A more general form of this expression is derived in Appendix 4, valid for any unbiased decoder.

Apart from producing an asymmetric $\mathrm{CP}(p_{\text{CR}})$ pattern, the gain fluctuations also create a negative covariation between the CP at $p_{\text{CR}}=0.5$ and the degree of asymmetry of the $\mathrm{CP}(p_{\text{CR}})$ pattern. This covariation appears because the cells with a higher portion of their variability driven by the gain fluctuations (higher ${\lambda }_i$) have a higher attenuation of ${\mathrm{CC}}_i(p_{\text{CR}}=0.5)$, given Equation 23a, while both a higher ${\lambda }_i$ and smaller ${\mathrm{CC}}_i(p_{\text{CR}}=0.5)$ lead to an increase in the slope ${\displaystyle {\beta }_{p_{\mathrm{C}\mathrm{R}}}\equiv {\sigma }_G{\lambda }_i(1-{\mathrm{C}\mathrm{C}}_i^2(p_{\mathrm{C}\mathrm{R}}=0.5))}$ of the dependence on $\mathrm{\Delta }s(p_{\text{CR}})$ in Equation 23b. Furthermore, a smaller ${\mathrm{CC}}_i(p_{\text{CR}}=0.5)$ also leads to a smaller effect of the multiplicative symmetric modulation $h(p_{\text{CR}})$, further contributing to the negative covariation between the magnitude of the CP and predominance of the symmetric or asymmetric pattern.

To illustrate the properties common to the model and to the $\mathrm{CP}(p_{\text{CR}})$ patterns from the MT data, Figure 4E shows CPs from Equation 23 as a function of $p_{\text{CR}}$ for examples combining four values of ${\mathrm{CC}}_{i0}(p_{\text{CR}}=0.5)$ and two values of ${\sigma }_G^2$, while the other parameters of the cell responses are kept constant. In particular, to determine ${\lambda }_i^2$ only in terms of the strength of the gain we fixed the rate to $f_i(p_{\text{CR}}=0.5)=10\mathrm{spike}/\mathrm{s}$ and considered the variance not associated with the gain to be equal to that rate, so that ${\lambda }_i^2=1/(1+0.1/{\sigma }_G^2)$. Accordingly, the values of ${\sigma }_G^2$ in Figure 4E, namely ${\sigma }_G^2=0.1$ and ${\sigma }_G^2=0.01$, correspond to ${\lambda }_i^2=0.5$ and ${\lambda }_i^2=0.09$, respectively. Further analysis of the model is provided in Appendix 4, where we also discuss the form of the $\mathrm{CP}(p_{\text{CR}})$ pattern produced by gain fluctuations when the decoder is composed by two pools of opposite choice preference (Shadlen et al., 1996).

To experimentally estimate the coefficients ${\beta }_{p_{\text{CR}}}^{(exp)}$ we fitted a quadratic regression of the CPs on the stimulus levels. To theoretically estimate the coefficients ${\beta }_{p_{\text{CR}}}^{(th)}$, we used the negative binomial model of Goris et al., 2014 to estimate ${\sigma }_G^2$ for each cell and used the form ${\displaystyle {\beta }_{p_{\mathrm{C}\mathrm{R}}}\equiv {\sigma }_G{\lambda }_i(1-{\mathrm{C}\mathrm{C}}_i^2(p_{\mathrm{C}\mathrm{R}}=0.5))}$ predicted by the gain model (Equation 11) to estimate ${\beta }_{p_{\text{CR}}}^{(th)}$.

We implemented a new GLM, called stimulus-dependent-choice GLM, that includes regression coefficients quantifying the effect on the firing rate of interactions between stimulus and choice. This model of the firing rate of each neuron was compared to two simpler and traditional models: a stimulus-only GLM, which includes only stimulus predictors of the neuron’s firing rate, and a stimulus-independent-choice GLM, which includes together with the stimulus predictor a single, stimulus-independent choice predictor.

In more detail, all three GLMs were Poisson models in which the mean firing rate $\mu (r_i)$ of cell $i$ was generally expressed by the following equation:

The terms ${\mathrm{\Sigma }}_{j=0}^4{a}_j{s}^j$ are present in all three types of GLM, and model the stimulus influence with a fourth order polynomial function of the stimulus level. These are the only terms of the stimulus-only GLM.

The choice dependence is modeled by ${\mathrm{\Sigma }}_{j=1}^{N_c}{I}_{P_j}(p_{\text{CR}})b_{j}D$, with the parameter $N_c$ ($N_c\in \{1,2,3\}$) setting the number of possible different levels of stimulus-dependent choice (we restricted the fitting to up to three different choice levels for simplicity, and because we found empirically this to work well for the MT data analyzed here). $I_{P_j}(p_{\text{CR}})$ is an indicator function which equals one if a $p_{\text{CR}}$ value belongs to the subset $P_j$ of values selected to be associated with the choice parameter b_j, and is zero otherwise. For the stimulus-independent-choice model, we set $N_c=1$ so that the choice affects the predicted responses equally for all stimulus levels. For the stimulus-dependent-choice GLM, we set $N_c>1$. For this stimulus-dependent-choice GLM, we determined the subsets of stimulus levels associated with each of those parameters using the $\mathrm{CP}(p_{\text{CR}})$ profiles for a first characterization of the stimulus dependencies. Like for the CP analysis, for each cell we determined which coherence values could be included in the analysis given a criterion requiring a minimum number of trials for each choice (at least 4). The existence of non-monotonic $\mathrm{CP}(p_{\text{CR}})$ profiles, such as the symmetric pattern around $p_{\text{CR}}=0.5$, indicated that it would be sub-optimal to tile the domain of $p_{\text{CR}}$ with $N_p$ bins and assign a different choice-parameter level to each bin. Accordingly, we first estimated the $\mathrm{CP}(p_{\text{CR}})$ profile of each cell and then used $k$-means clustering with an Euclidean distance to cluster the components of $\mathrm{CP}(p_{\text{CR}})$, corresponding to the bins of $p_{\text{CR}}$, into $N_c$ subsets. A different GLM choice-parameter b_j was then assigned to each choice-parameter level $j=1,\mathrm{\dots }{N}_c$.

We compared the predictive power of the three types of models using cross-validation. To avoid that the choice-parameters fitted were affected by the ratio of trials with each choice, we matched the number of trials of each choice used to fit the model at each choice-parameter level. We first merged in two pools, one for each choice separately, the trials of all stimulus levels assigned to the same choice-parameter level. We then determined the number of trials from each pool to be included in the fitting set as an 80% of the trials available in the smallest pool, hence matching the number of trials selected from each choice. The remaining trials were left for the testing set. This procedure was repeated for each choice-parameter level and a GLM model was fitted on the fitting set obtained combining the selected trials for all levels. This random separation between fitting and testing data sets was repeated 50 times and the average predictive power was calculated. Performance was then quantified comparing the increase in the likelihood of the data in the testing set with respect to the likelihood of the null model which assumes a constant firing rate (L₀). To determine if incorporating the choice as a predictor improved the prediction, we examined the relative increase in likelihood (RIL) defined as the ratio of the likelihood increase $L(choice,stimulus)-L(stimulus)$ and the increase $L(stimulus)-L_0$. For the stimulus-dependent-choice models, we selected the most predictive model from $N_c=2,3$. To evaluate the improvement when considering stimulus-dependent choice influences, we compared the RIL obtained for the stimulus-dependent-choice and stimulus-independent-choice models.

The codes for the analysis of Choice Probability stimulus-dependencies and GLMs with stimulus-choice interaction terms are available at https://github.com/DanielChicharro/CP_DP (copy archived at swh:1:rev:5850c573860eb04317e7dc550f96b1f47ca91c6a).

We here provide details of how the CP analytical expression of Equation 16 is obtained from the definition of Choice Probability (Equation 12) when the probability of the responses for each choice has the form of Equation 15, derived from the threshold model. We subsequently characterize the statistical power for the detection of the threshold-induced modulation $h(p_{\text{CR}})$, as a function of the magnitude of the CP, the number of trials, and the number of cells used to estimate an average $\mathrm{CP}(p_{\text{CR}})$ profile.