This section provides an overview of the Quadratic Color Model (QCM); a full mathematical description is provided in the Appendix 1. Given a description of the stimulus, the QCM provides a prediction of the BOLD fMRI response within V1. Our stimuli were full field temporal chromatic modulations that can be specified by their contrast (vector length of the stimulus in the LM contrast plane) and chromatic direction (angle of the stimulus measured counterclockwise with respect to the positive abscissa). From this stimulus specification, the model employs three stages that convert the input to the BOLD fMRI response (Figure 1). First, a quadratic isoresponse contour is defined that allows for the transformation of contrast and direction into what we term the ‘equivalent contrast’. Second, a single non-linear function transforms the equivalent contrast to a prediction of the population neuronal response underlying the BOLD response. Finally, the neuronal response is converted to a predicted BOLD response by convolution with a hemodynamic response function.

Figure 1

Download asset Open asset

The first stage of the QCM computes the equivalent contrast of a stimulus from its cone contrast using a subject-specific elliptical isoresponse contour. Equivalent contrast is the effective contrast of a stimulus in V1 once it has been adjusted to account for differences in the neuronal sensitivity to stimulation across different chromatic directions. An isoresponse contour is defined as a set of stimuli that evoke the same neuronal response. In the QCM, the loci of such stimuli form an elliptical isoresponse contour in the LM cone contrast plane. All points on this elliptical isoresponse contour have the same equivalent contrast (Figure 1A, dashed gray ellipses). As the amplitude of the neuronal response increases, the ellipse that defines the set of stimuli producing that response also grows in overall scale. Importantly, the QCM assumes that the aspect ratio and orientation of elliptical isoresponse contours do not change as a function of the response level; only the overall scale of the ellipse changes. The use of elliptical isoresponse contours is motivated by prior psychophysical (Poirson et al., 1990; Knoblauch and Maloney, 1996), electrophysiological (Horwitz and Hass, 2012), and fMRI experiments (Engel et al., 1997) which have successfully used ellipses to model chromatic isoresponse contours.

The elliptical isoresponse contours are described by a symmetric quadratic function that defines the major and minor axes of the ellipse. We use this quadratic function to compute the equivalent contrast for each stimulus. The vector lengths of all stimuli that lie on a single isoresponse contour provide the cone contrasts required to elicit an equal neuronal response. The minor axis of the elliptical isoresponse contour corresponds to the chromatic direction that requires the least amount of cone contrast to produce this equal neuronal response, and is therefore the most sensitive chromatic direction. The major axis corresponds to the direction of least sensitivity. At this stage, the model is only concerned with the shape of the elliptical contour, thus we adopt the convention of normalizing the ellipse used to define equivalent contrast so that its major axis has unit length. This allows the length of the minor axis to directly represent the relative sensitivity, which is taken as a ratio of the minor axis (maximal sensitivity) to major axis (minimal sensitivity), referred to as the minor axis ratio. The angle of the major axis in the LM contrast plane (ellipse angle) orients these maximally and minimally sensitive directions.

Since all of the stimuli that lie on a single isoresponse contour produce the same response, we can represent these points by their common equivalent contrast. The neuronal responses to stimuli across different color directions are a function of this single variable, and therefore we can transform equivalent contrast into predicted neuronal response via a single static non-linear function (Figure 1B). Here, we employ the four-parameter Naka-Rushton function (see Appendix 1).

To predict the BOLD fMRI signal, we obtain the time-varying neuronal response prediction from the Naka-Rushton function for a stimulus sequence presented in the fMRI experiment. This neuronal response is convolved with a subject-specific hemodynamic response function to produce a prediction of the BOLD fMRI signal (Figure 1C).

In summary, the QCM takes as input the temporal sequence of stimulus modulations, defined by their chromatic direction and contrast in the LM cone contrast plane, and outputs a prediction of the BOLD fMRI time course. The QCM has six free parameters: two that define the shape of the normalized elliptical isoresponse contour and four that define the Naka-Rushton equivalent contrast response function.

We examined how well the GLM fit the measured BOLD response from area V1. Figure 3 shows the fit of the GLM for six example runs from Subject 2. In each panel, the measured BOLD percent signal change is shown as the thin gray line, while the fit obtained from the GLM is shown as the orange line. The orange shaded region represents the 68% confidence interval of the fit found using bootstrap resampling. The GLM fit captured meaningful stimulus-driven variation in the BOLD response, with some variation in fit quality across runs. The median R² value across runs was 0.41 for Measurement Set 1 and 0.32 for Measurement Set 2. Fits for the other two subjects are provided as Figure 3—figure supplements 1–2. Due to the randomized stimulus order within each run, it was not straightforward to determine the degree to which the unmodeled variance was due to stimulus-driven structure not modeled by the GLM (e.g. carry-over effects) as opposed to measurement noise. Overall, the quality of the GLM fits supported using the GLM as a benchmark model, as well as using the GLM beta weights as a measure of the V1 response.

Figure 3 with 3 supplements see all

Download asset Open asset

We used a leave-runs-out cross-validation procedure to compare the GLM and the QCM (see Materials and methods section for details). This cross-validation compares the ability of the models to predict data not used to fit the parameters, accounting for the possibility that more flexible models (such as the GLM) may overfit the data. Figure 4 shows the results of the cross-validation comparison for all subjects. Both models track meaningful variation in the signal, although less so for the data from Subject three in Measurement Set 2. Importantly, we see that the QCM cross-validated R² is essentially indistinguishable from the GLM cross-validated R², although in all cases slightly higher.

Figure 4 with 1 supplement see all

Download asset Open asset

To further assess differences between the GLM and the QCM, we analyzed the model residuals as a function of the stimulus condition, to check for systemic patterns in the residuals as well as any differences between the two models in this regard. We first examined each direction/contrast pair separately by plotting the residuals of the GLM and QCM over the 14 TRs after the start of each stimulus block. We did not observe any systematic variation in the residuals as a function of contrast level within a single chromatic direction. Therefore, we examined the mean residual value, taken from 4 to 14 TRs after stimulus onset, for all trials in a chromatic direction (collapsed over contrast). We plot these mean residuals for both the GLM and the QCM, as a function of chromatic direction, for each subject and session in Figure 3—figure supplement 3. From this, we observe no consistent pattern of residuals within or across models. Note that the residual values for the GLM and QCM mostly overlap, despite the GLM having separate parameters for each stimulus direction.

We also employed a leave-session-out cross-validation procedure to assess the generalizability of the QCM (See Materials and methods for details). Given that Sessions 1 and 2 do not share any common chromatic directions, we were able to evaluate how effectively the QCM generalizes to chromatic directions not used to derive the model parameters. The green contrast response functions shown in Figure 5 result from fitting the QCM to either Session 1 or Session 2, and predicting the responses from the held-out session. The generalization from Session 1 to Session 2 (right-hand subplots) is excellent for this subject. The generalization from Session 2 to Session 1 is also good, albeit with a large confidence interval for the 45° direction. For other subjects and measurement sets, the QCM generalizes reasonably well (Figure 5—figure supplements 1–3). Overall, generalizations from Session 1 to Session 2 perform better than those from Session 2 to Session 1. This finding may reflect the particular set of chromatic directions presented in each session: only Session 1 includes a chromatic direction close to the major axis of the ellipse, which better constrains the QCM fit. Therefore, the QCM is capable of generalizing well to unmeasured chromatic directions, with the requirement that the stimuli include chromatic directions and contrasts that adequately constrain the model parameters.

Figure 5 with 3 supplements see all

Download asset Open asset

Conceptually, the parameters of the QCM characterize two key model components. The first component defines the contrast-independent shape of elliptical isoresponse contour. This describes the relative sensitivity of V1 to modulations in all chromatic directions within the LM contrast plane. The second component defines the response nonlinearity, which is independent of chromatic direction. It operates on equivalent contrast to produce the underlying neural response.

The isoresponse contour is described by two parameters: the direction of least sensitivity (ellipse angle; counterclockwise to the positive abscissa) and the ratio of vector lengths between the most and least sensitive directions (minor axis ratio; see Quadratic Color Model Section and Appendix 1). Within the QCM, the angle and minor axis ratio provide a complete description of chromatic sensitivity that is contrast independent.

Figure 6 shows the QCM isoresponse contours for all three subjects and both measurement sets. We found that for all subjects and measurement sets, the angle of the isoresponse contours was oriented at approximately 45°. An ellipse angle of 45° indicates that V1 was least sensitive to stimuli modulated in the L+M direction, and most sensitive to stimuli modulated in the L-M direction. Across all subjects and measurement sets, the minor axis ratio parameters ranged between 0.15 and 0.25. Thus, for the spatial and temporal properties of our modulations, V1 was roughly five times more sensitive to modulations in the L-M direction than the L+M direction. We found good agreement between the isoresponse contours from the independent measurement sets as well as across subjects. 

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 7 shows the V1 equivalent contrast nonlinearity of the QCM for Subject 2 for both measurement sets. This non-linearity describes how the underlying neuronal response increases with increasing equivalent contrast. We used the isoresponse contour of the QCM to convert the chromatic direction and cone contrast of each stimulus to its equivalent contrast. This allowed us to replot each beta weight derived from the GLM (Figure 2) on an equivalent contrast axis (Figure 7; closed circles). For all subjects, the single nonlinearity accurately captured the dependence of the GLM beta weight on equivalent contrast, with no apparent bias across chromatic directions. The agreement between the GLM beta weight points and QCM fits demonstrated that separating the effects of chromatic direction and contrast in the QCM is reasonable. Figure 7—figure supplement 1 provides the same plots for Subjects 1 and 3. 

Figure 7 with 1 supplement see all

Download asset Open asset

We leveraged the QCM to examine how chromatic sensitivity varies with eccentricity within V1. Figure 9 plots the V1 QCM parameters as a function of eccentricity, for Subject 2. The left panel shows the minor axis ratio and the right panel shows the ellipse angle. In both plots, individual points represent a single vertex, with the x-axis giving the visual field eccentricity of that vertex obtained from the Benson et al., 2014 template, and the y-axis giving the parameter value. The transparency of each point indicates the R² value of the QCM fit for the corresponding vertex. The maximum R² value across vertices for Measurement Sets 1 and 2 were 0.25 and 0.24, respectively. The lines in each panel reflect a robust regression fit to the points. We found that there is little change in either parameter with eccentricity. For Subject 2, the best fit lines had slightly negative slopes for the minor axis ratio and slightly positive slopes for the ellipse angle, with good agreement across measurement sets. The overall change in parameter values from 0° to 20°, however, was small compared to the vertical spread of values at each eccentricity. We compared the change in parameter values from 0° to 20° to the variability across measurement sets for all three subjects (Table 1 for minor axis ratio, and Table 2 for ellipse angle). Across subjects, the majority of sessions showed small differences in parameter values from 0° to 20°, but we note that these did in some cases exceed the measurement set-to-set difference in the parameter values obtained for all of V1.

Figure 9 with 2 supplements see all

Download asset Open asset

The plots shown in Figure 9 examine how the QCM parameters vary with eccentricity. To allow comparison with prior studies of how the BOLD response varies with eccentricity within V1 (Vanni et al., 2006; Mullen et al., 2007; D'Souza et al., 2016), we also used the QCM to predict how the response would vary for stimuli in the L+M and the L-M directions, and plot these predicted responses as a function of eccentricity (Figure 9—figure supplement 1). This was done on a vertex-by-vertex basis, using within-vertex QCM parameters. Specifically, we chose the 50% contrast stimulus condition for both the L-M and the L+M direction (contrasts of 0.06 and 0.30 respectively). Using these stimuli, we computed the predicted neuronal response by applying the QCM forward model (the transformation to equivalent contrast and the Naka-Rushton steps) using the parameter values corresponding to that particular vertex. Examining the data in this way reveals a negligible change in response as a function of eccentricity for both the L+M and L-M directions, for all subjects and measurement sets. We return in the discussion to consider the relation between our results and those found in prior studies.

A goal of systems neuroscience is to link measurements of neuronal properties to measurements of behavior. To make these links, the measurements made in each domain must be placed into a common space for comparison. The QCM provides a way to represent fMRI measurements in a manner that makes such comparisons straightforward. The contrast-invariant isoresponse contour from the QCM provides us with a stimulus-referred characterization of the BOLD fMRI response. Other methodologies, such as psychophysics or electrophysiology, may be used to obtain similar characterizations, allowing for comparisons across response measures within this common framework. For example, an approach to studying chromatic sensitivity is to characterize the isothreshold contour, which specifies the set of stimulus modulations that are equally detectable. Engel et al., 1997 took this approach and found that for low temporal frequencies the psychophysical isothreshold and BOLD fMRI isoresponse contours in the LM contrast plane were well-described as ellipses and had similar shapes. While some work has argued that psychophysical isothreshold may deviate subtly from ellipses (for review see Stockman and Brainard, 2010), two studies that attempted to reject the elliptical form of such contours did not do so (Poirson et al., 1990; Knoblauch and Maloney, 1996). Consistent with Engel, Zhang, and Wandell (1997), we found elliptical BOLD isoresponse contours at our 12 Hz temporal frequency with highest sensitivity in the L-M direction. As they note, although psychophysical isothreshold contours remain well-described by ellipses at higher temporal frequencies, their orientation changes to favor L+M sensitivity over L-M sensitivity. This dissociation in the particulars of the isothreshold and BOLD isoresponse contours makes it unlikely that the mechanisms that contribute to the BOLD response in V1 limit psychophysical detection at the higher temporal frequencies, unless there are important temporal-frequency dependent changes in response variability that are not captured by the BOLD measurements.

Many theories of color vision postulate that signals from the L-, M-, and S-cone photoreceptors are combined to form three post-receptoral mechanisms, roughly characterized as an additive combination of L- and M-cone contrast (L+M), an opponent combination of L- and M-cone contrast (L-M), and an opponent combination of S-cone contrast with L- and M-cone contrasts (S-(L+M)) (Stockman and Brainard, 2010; Shevell and Martin, 2017). Our finding that the major and minor axes of the isoresponse ellipse are well-aligned with the L+M and L-M modulation directions agrees with such theories. More generally, a quadratic isorepsonse contour can be produced by a quadratic mechanism that computes a sum of the squared responses of two underlying linear mechanisms, where the output of each linear mechanism is a weighted sum of L- and M-cone contrasts (Poirson et al., 1990). If the two linear mechanisms are L+M and L-M mechanisms with the weights appropriately chosen to represent the relative sensitivities (L-M sensitivity greater than L+M sensitivity), then the isoresponse contour of the resulting quadratic mechanism will be a close match to those we measured.

Note, however, that other pairs of underlying mechanisms are also consistent with the same elliptical isoresponse contours (Poirson et al., 1990), so that our isoresponse contours do not uniquely determine the sensitivity of the underlying linear mechanisms, even within the QCM together with the assumption that there are two such mechanisms.

More generally, one can construct non-quadratic models whose isoresponse contours approximate the ellipse we found using the QCM, and if this approximation is good our data will not reject such models. To illustrate this point, we developed and fit an alternate model, the Linear Channels Model (LCM), a variation on the Brouwer and Heeger channel model (Brouwer and Heeger, 2009; Kim et al., 2020), that accounts for our data about as well as the QCM (see Appendix 1; Figure 4—figure supplement 1). The best fitting isoresponse contours found with the LCM, which could in principle deviate considerably from an ellipse, none-the-less approximate the ellipse we found using the QCM, but are not perfectly elliptical (Figure 6—figure supplement 1). Despite the agreement at the functional level of the isoresponse contours, the properties of the mechanisms underlying the LCM differ from those of the QCM, and these properties also differ across different instantiations of the LCM that account for the data equally well (see Appendix 1).

Because of the similarity in cross-validated R² values across the GLM, Naka-Rushton, QCM and LCM models, the reader may wonder whether the data have sufficient power to reject any isoresponse contour shape. To address this, we fit and cross-validated a form of the QCM with the angle constrained to 0 degrees. This resulted in a noticeably lower cross-validated R² for this model as compared to all other models we tested (Figure 4—figure supplement 1, labeled at ‘QCM locked’) and provides reassurance that the data indeed have power to inform as to the shape of the isoresponse contour. We expect that other isoresponse contour shapes that differ from the best fitting QCM contour to a degree similar to that of the constrained ellipse would also be rejected.

Thus, while measurements of the BOLD response place constraints on the population response properties of the neuronal mechanisms, these properties are not uniquely determined given the BOLD response alone. The ambiguity is further increased if we consider properties of individual neurons, as the aggregate BOLD response will be shaped both by the response properties of such neurons and the numbers of different types of neurons in the overall neural population. With that caveat, we make some general observations. Our stimuli were large, spatially uniform (effectively 0 c.p.d.) chromatic modulations that were temporally modulated at 12 Hz. These stimuli could drive ‘color’, ‘luminance’ and ‘color-luminance’ cells as described by Johnson et al., 2001, depending on the particular chromatic direction. The 0 c.p.d. stimuli would presumably produce strong responses in the ‘color’ cells for our L-M direction, given that these cells are thought to behave as low-pass filters in the spatial domain and have unoriented receptive fields. It is less clear how strongly ‘luminance’ and ‘color-luminance’ cells would respond to our spatially uniform stimuli given that these cells are spatially bandpass. Schluppeck and Engel, 2002 plot the spatial frequency response function estimated from the data from Johnson et al., 2001. These functions plot the average firing rate as a function of spatial frequency for the color, color-luminance, and luminance cells. Taking the lowest spatial frequency present in the dataset (0.1 c.p.d.), firing rates for color cells are roughly 5x times higher than the firing rates for luminance and color-luminance cells. This is the same as our average minor axis ratio which indicates V1 is roughly five times more sensitive to L-M than to L+M. A caveat here is that the 12 Hz flicker rate of our stimuli might shift responses relative to the analysis of Schluppeck and Engel, 2002.

It is also important to note that our data do not distinguish the extent to which the response properties of the BOLD signals we measure in V1 are inherited from the LGN or are shaped by processing within V1. The spatiotemporal properties of our stimuli would robustly drive cells in cell in the LGN that project to V1 (Lankheet et al., 1998). Even though the signals measured in our experiment are spatially localized to V1, we cannot ascribe the observed response properties to particular neural processing sites. As such, the V1 sensitivities found from fitting the QCM may be inherited from areas prior to V1. In principle, they could also be affected by feedback from other cortical areas. We also undertook an analysis of data from the LGN, but found that these signals were too noisy to reveal reliable stimulus-driven responses in our data.

Finally, it is interesting to observe that quadratic models have been used to characterize the isoresponse properties of individual neurons in macaque V1 (Horwitz and Hass, 2012). Roughly half of the neurons tested in this paper were best fit by quadratic isoresponse surfaces (in the L-, M-, and S-cone contrast space) while the other half were well fit by a linear model whose isoresponse surfaces were parallel planes. Of the quadratic isoresponse surfaces, some were ellipsoidal while others were hyperbolic. As noted above, despite this qualitative similarity, connecting the diverse population of individual neural responses to the aggregated BOLD fMRI response remains a challenge for future work. The QCM fit to our data aids in this endeavor through the constraint its isoresponse contour places on the aggregated neural response.

Many aspects of visual function change with eccentricity (Rosenholtz, 2016), and understanding and quantifying this variation is a key part of a functional characterization of vision. In addition, prior work attempts to relate such functional variation with eccentricity to variation in the underlying neural mechanisms. Relevant to the present work is the idea that variation of chromatic sensitivity with eccentricity can inform as to how signals from separate cone classes are combined by retinal and cortical neural circuitry (Lennie et al., 1991; Mullen and Kingdom, 1996; Wool et al., 2018; Baseler and Sutter, 1997). In this context, we examined how the parameters of the QCM for individual vertices varied with eccentricity across V1. Overall, we find little change in the isoresponse contours with eccentricity within the central 20° of the visual field (Figure 9; Tables 1 and 2; Figure 9—figure supplement 2), and an analysis of how predicted response to L+M and L-M modulations would change with eccentricity also shows little or no effect (Figure 9—figure supplement 1). In addition, we do not observe any clear change in the parameters of the contrast-response function with eccentricity (analysis not shown). Overall, the BOLD response to chromatic modulations, as evaluated in our data using the QCM, is remarkably stable across V1. That the orientation of the elliptical isoresponse contours does not change with eccentricity is consistent with psychophysical studies that show that the relative L- and M-cone inputs to an L-M mechanism are stable across the visual field (Newton and Eskew, 2003; Sakurai and Mullen, 2006). The fact that our data do not show a loss in L-M sensitivity relative L+M sensitivity with increasing eccentricity, on the other hand, is not commensurate with psychophysical studies that do show such a loss (Stromeyer et al., 1992; Mullen and Kingdom, 2002; Mullen et al., 2005; Hansen et al., 2009). As discussed above, BOLD fMRI sensitivity in V1 does not always mirror psychophysical sensitivity. Thus we focus below on comparison between our results and other fMRI studies of how sensitivity in V1 varies with eccentricity.

Several prior studies have used BOLD fMRI to examine how visual cortex responses vary with eccentricity to stimuli modulated in L-M and L+M directions (Vanni et al., 2006; Mullen et al., 2007; D'Souza et al., 2016). Both Mullen et al., 2007 and Vanni et al., 2006 report a decrease in the V1 response to L-M modulations with eccentricity, while the response to L+M remains roughly constant. This differs from our result, and the size of the effects in these papers are large enough that we would expect that if they were present in our data they would be visible in the analysis shown in Figure 9—figure supplement 1 (see Figure 8 in Mullen et al., 2007 and Figure 8 in Vanni et al., 2006). Both speculate that the mechanism underlying this observation is non-selective (random) connections between L and M cones and retinal ganglion cell receptive fields. If these connections are non-selective, L-M sensitivity would be expected to decrease with eccentricity. This is because the area in which receptive fields pool cone inputs increases with distance from the fovea, progressively reducing the likelihood that random L- and M-cone inputs to the center and surround will produce chromatic opponency (Lennie et al., 1991; Mullen and Kingdom, 1996; Wool et al., 2018). In contrast, D'Souza et al., 2016 find, for the majority of spatial frequencies studied, no change in L-M response relative to the response to isochromatic luminance modulations. This result is generally in line with our data, although at their lowest spatial frequencies D’Souza and colleagues observe a modest decline in relative L-M sensitivity. Following the same line of reasoning as Mullen et al., 2007 and Vanni et al., 2006 but reaching the opposite conclusion, D'Souza et al., 2016 take their result as supporting the idea that connections between cones and some classes of ganglion cells are selective for cone type and preserve chromatic sensitivity across the retina.

Comparison across our and the prior studies is complicated by variation in the stimuli used. Indeed, a dependency on spatial frequency is indicated by the data of D'Souza et al., 2016. We used spatially uniform fields, while other studies use stimuli with higher-spatial frequency content. More generally, other factors could also lead to variation across studies, as well as complicate inferring the properties of retinal wiring from how psychophysical thresholds or measurements of cortical response vary with eccentricity.

One such factor is the changes in cone spectral sensitivity with eccentricity, caused primarily by variation in macular pigment and photopigment optical density. Variation in macular pigment and photopigment optical density can produce eccentricity-dependent deviations in the degree of actual cone contrast reaching the photoreceptors. Prior studies do not account for this variation, leading to the possibility that effects of eccentricity on sensitivity are due to receptoral, rather than post-receptoral mechanisms. In our study, we designed spectral modulations that produce the same contrasts in both 2-degree and 15-degree cone fundamentals (see Materials and methods and Tables 5–7). This reduces the change in cone contrast with eccentricity for our stimuli.

Another factor, not emphasized in previous work, is that the size of the effects will depend on where on the underlying contrast response functions the responses to the stimuli in the chromatic directions being compared lie. To understand this issue, consider the example in which the response for one chromatic direction is well into the saturated regime of the contrast-response function, while the other direction is not. This could lead to artifacts in the measured ratio of the response to the two directions with eccentricity, where any change in response for the saturated direction is hidden by a ceiling effect. Our study minimizes the role of this factor through measurement and modeling of the contrast response functions in each chromatic direction, allowing the QCM to extract a contrast independent shape for the isoresponse contours.

Although changes in sensitivity with eccentricity can be caused by mechanisms at many levels of the visual system, the lack of such variation in our data is parsimoniously explained by retinal output that preserves chromatic sensitivity with eccentricity. This interpretation is challenged by studies that show random (Wool et al., 2018) or close to random (Field et al., 2010) inputs from L- and M-cones to midget ganglion cell centers in the retinal periphery. Not all studies of midget cell chromatic responses or their parvocelluar LGN counterparts agree with non-selective wiring (Reid and Shapley, 1992; Martin et al., 2001; Martin et al., 2011; Lee et al., 2012). One possible cortical mechanism that could compensate for the reduced L-M signal-to-noise ratio with eccentricity as predicted by random wiring models is supra-threshold compensation for reduced signals. Mechanisms of this sort have been postulated in the domain of contrast perception (Georgeson and Sullivan, 1975) and anomalous trichromatic vision (Boehm et al., 2014; Tregillus et al., 2021). Another possibility is a differential change in stimulus integration area across chromatic directions and with eccentricity. However, this latter possibility is not supported by fMRI population receptive field (pRF) measurements for modulations in different chromatic directions in V1 (Welbourne et al., 2018).

In our study, we only measured responses to stimuli confined to the LM contrast plane. A more general account of chromatic contrast sensitivity requires modulating stimuli in all three-dimension of the full L-, M-, and S-cone contrast space. The QCM may be generalized in a straightforward manner to handle this expanded stimulus set by replacing the elliptical isoresponse contours with ellipsoidal isoresponse surfaces, but we have yet to test this generalization.

Another way in which the QCM may be generalized, even within the LM contrast plane, is to consider modulations at other temporal frequencies. In our experiment, we fixed the temporal modulation of the stimulus at 12 Hz. The QCM could be readily fit to data from modulations at various other temporal frequencies with the goal of observing how the chromatic sensitivity changes. With this, one could further examine the BOLD response to mixtures of different temporal frequency modulations. This is particularly interesting in that any arbitrary complex temporal modulation can be decomposed into an additive mixture of modulations at different temporal frequencies and phases (Bracewell, 1978). If the response to temporal frequency mixtures can be predicted via a simple rule of combination (such as linearity), establishing the QCM parameters for a well-chosen set of temporal frequencies would enable prediction of the BOLD response to chromatic stimuli modulated with complex temporal sequences.

Just as we had a fixed temporal frequency in our experiment, the stimulus presented also had a fixed spatial frequency (0 c.p.d.). Similar to the temporal domain, complex spatial images can be broken down into a combination of oriented two-dimensional sine wave patterns at single spatial frequencies and phases. Examining how the QCM fits change with changing spatial frequencies might allow for models of the BOLD response to arbitrary complex spatial stimuli. Consistent with this general goal, recent work has developed quantitative forward models of the BOLD response of early visual cortex to a variety of achromatic modulations with different spatial patterns (Kay et al., 2013a; Kay et al., 2013b). These models have sequential stages of processing that operate on an input image and transform it into a model of the BOLD response. How such models should be generalized to handle chromatic modulations is not known. If the QCM holds for other spatial frequencies, such a result would place important constraints on the appropriate generalization for such forward models to incorporate color.

Three subjects (age 23, 25, and 26 years; two female) took part in the fMRI experiment. All subjects had normal or corrected to normal acuity and normal color vision. The research was approved by the University of Pennsylvania Institutional Review Board. All subjects gave informed written consent and were financially compensated for their participation.

Each subject participated in four sessions of data collection. The first two sessions constituted Measurement Set 1, and the second two sessions Measurement Set 2 (Measurement Set 2 being a replication of Measurement Set 1). In both sessions, subjects underwent 48 min of fMRI scanning, with Session 1 also including two anatomical scans (a T1-weighted and a T2-weighted scan). Subjects were tested for color vision deficiencies in a separate session, using the Ishihara pseudoischromatic plates (Ishihara, 1977). All subjects passed with no errors. The experimental procedures for Measurement Set 1 were preregistered (https://osf.io/wgfzy/), and an addendum describes the replication Measurement Set 2 (https://osf.io/zw6jp/).

All stimuli were generated using a digital light synthesis device (OneLight Spectra). This device produces desired spectra through the use of a digital micro-mirror device chip that in the configuration we used allows for the mixture of 56 independent primaries with a FWHM of ~16 nm and a refresh rate of 100 Hz.

Stimuli were generated to evoke specific photoreceptor responses through the use of silent substitution (Estévez and Spekreijse, 1982). Silent substitution operates on the principle that there exist sets of light spectra that, when exchanged, selectively modulate the activity of specified cone photoreceptors. Thus, stimulus modulations relative to a background can be generated such that they nominally modulate the activity of only the L-, M-, or S-cones, or combinations of cone classes at specified contrasts. Additional information on how the stimuli were generated is provided in Spitschan et al., 2015.

The stimuli account for differences in the cone fundamentals between the fovea and the periphery. This was done by treating the L, M, and S cones in the 2- and 15-degree CIE physiologically-based cone fundamentals (CIE, 2007) as six classes of photoreceptors. For any desired set of L-, M-, and S-cone contrasts, we designed modulations that attempted to provide the same contrasts on the 2- and 15-degree L-cone fundamentals, on the 2- and 15-degree M-cone fundamentals, and on the 2- and 15-degree S-cone fundamentals. This is possible because our device has 56 primaries, rather than the typical 3 of RGB displays. Our procedure has the effect of creating light spectra that reduce differences in the L-, M-, and S-cone contrasts produced across the retina. The cone fundamentals were tailored to the age of each subject, to account for age-related differences in typical lens density (CIE, 2007). See Tables 5–7 for the central and peripheral maximum stimulus contrast values for each subject and measurement set.

Spectroradiometric measurements of the stimuli were made before and after each experimental session. During the measurements made prior to the experiment, a correction procedure was run in which the spectral power distribution of the modulation in each chromatic direction were adjusted to minimize the difference between the measured and desired cone contrasts for the 2- and 15-degree cone fundamentals. These corrections were made to the modulation spectra at the maximum contrast used in each direction. The order in which spectroradiometric measurements were taken during an experimental session was (1) five pre-correction measurements for each chromatic direction used in the session, (2) the corrections procedure, (3) five post-correction measurements per direction, and (4) five post-experiment measurements per direction. The mean of post-correction and post-experiment cone contrast measurements for the individual subjects and measurement sets are provided in Tables 5–7.

The stimuli were confined to the LM plane of cone contrast space (Figure 10A; see also Figure 1A). Cone contrast space has three axes that are defined by the relative change in the quantal catch of the L, M, and S cones when modulating between the light spectra of interest and a specified reference spectrum. We refer to the reference spectrum used to calculate this relative change in cone excitations as the background (nominal chromaticity; x = 0.459, y = 0.478, luminance Y = 426 cd/m^2; chromaticity and luminance computed with respect to the XYZ 10° physiologically-relevant color matching functions, as provided at https://cvrl.org). The background corresponds to the origin of cone contrast space. The LM contrast plane is a subspace of cone contrast space consisting of modulations that affect only L- and M-cone excitations, but which leave S-cone excitations unchanged relative to the background. A point in cone contrast space specifies how much L- and M-cone contrast is produced by modulating from the background to the specified stimulus. Points lying along the x-axis of Figure 10A modulate only L-cone contrast while M-cone contrast remains constant. Points lying along the y-axis modulate only M-cone contrast while keeping L-cone contrast constant. Points in intermediate directions modulate both L- and M-cone contrast, in proportion to the x- and y-axis components.

Figure 10

Download asset Open asset

All stimuli were spatially uniform (0 c.p.d.), full-field temporal modulations with a radius of 30° visual angle. The temporal modulations were bipolar sinusoids around the reference background, with the positive and negative arms designed to increase and decrease targeted cone excitations in a symmetric fashion. All stimuli were modulated at 12 Hz (Figure 10B). A single modulation is thus described by pair of vectors in the LM contrast plane that have an angle of 180° between them, corresponding to the positive and negative arms of the modulation (Figure 10A). The entries of the vectors are the L- and M-cone contrasts of the end points of each arm. In this paper, we refer to a modulation by the angle made between its positive arm and the positive x-axis (corresponding to 0° in the LM contrast plane), with angle increasing counterclockwise. We refer to each angle tested as a chromatic direction. In total, we tested the eight chromatic directions: −45°, −22.5°, 0°, 22.5°,45°, 67.5°, 90°, and 112.5°. The −45°, 0°, 45°, and 90° directions correspond to L-M, L-cone isolating, L+M, and M-cone isolating directions, respectively. For all chromatic directions, the spectra were designed to produce constant S-cone contrast.

We express the stimulus contrast of a modulation as the vector length (L2 norm) of the cone contrast representation of its positive arm (Figure 10A). Gamut limitations of the light synthesis engine result in different maximum stimulus contrasts for different chromatic directions (see heavy dashed contour in Figure 10A). The maximum contrast used in each direction is provided in Table 3. For all directions, we tested five contrast levels. The contrast levels tested for each chromatic direction were selected to be log spaced relative to the maximum contrast used. The relative contrasts were 100, 50, 25, 12.5, and 6.25 percent. We also measured a 0 contrast reference condition in which the background without modulation was presented.

We measured whole brain BOLD fMRI responses to stimuli modulated in eight different chromatic directions, each with five contrast levels, using a block design. In total, we tested 40 different combinations of contrast levels and chromatic directions. We split the eight chromatic directions into two separate scanning sessions of four directions each (Figure 11A). In Session 1, we tested −45° (L-M), 0° (L Isolating), 45° (L+M), and 90° (M Isolating). In Session 2, we tested the other four directions (−22.5°, 22.5°, 67.5°, and 112.5°). The order of data collection for the two sessions was randomized across subjects. A session consisted of 10 runs and each run had a duration of 288 s. Within a run, each contrast/direction pair was presented within 12 s blocks (Figure 11). The order of contrast/direction pairs was psuedorandomized within each run. Along with four presentations of a background-only block, each run consisted of 24 blocks. The background-only blocks contained no temporal contrast modulation, providing a reference condition for data analysis. We chose 12 Hz modulations based upon prior work showing that for stimuli similar to the ones used in this study (L+M and L-M), this frequency elicited a robust response in V1 (Spitschan et al., 2016). Within a block, modulations were ramped on and off using a 500 ms half-cosine. Figure 10C provides a schematic of the structure of an functional run.

Figure 11

Download asset Open asset

During each functional run, subjects engaged in an attention task. This task consisted of pressing a button every time the stimulus dimmed (Figure 10C). Each attentional event lasted for 500 ms. The probability of an attentional event occurring in a block was 66% in Measurement Set 1 and a 33% in Measurement Set 2. The onset time of an attentional event within a block was random except that the event could not occur during the on and off half-cosine ramp. The purpose of the attention task was to encourage and monitor subject wakefulness throughout the scan session. All subjects responded to 100% of the attentional events throughout all runs and sessions.

Each subject was studied during an initial pair of scanning sessions that we call Measurement Set 1 (Figure 11A), and a subsequent pair of replication scans that we call Measurement Set 2 (Figure 11B). Measurement Set 2 tested stimuli with the same chromatic directions and contrast levels as Measurement Set 1. The grouping of chromatic directions within a session was the same across measurement sets. The two measurement sets used different pseudo-randomized presentation orders. Both measurement sets also randomly assigned which session was acquired first. Across both sessions and measurements sets, we collected a total of 960 blocks per subject. The two measurement sets were analyzed separately.

MRI scans made use of the Human Connectome Project LifeSpan protocol (VD13D) implemented on a 3-Tesla Siemens Prisma with a 64-channel Siemens head coil. A T1-weighted, 3D, magnetization-prepared rapid gradient-echo (MPRAGE) anatomical image was acquired for each subject in axial orientation with 0.8 mm isotropic voxels, repetition time (TR) = 2.4 s, echo time (TE) = 2.22 ms, inversion time (TI) = 1000 ms, field of view (FoV) = 256 mm, flip angle = 8°. BOLD fMRI data were obtained over 72 axial slices with 2 mm isotropic voxels with multi-band = 8, TR = 800 ms, TE = 37 ms, FOV = 208 mm, flip angle = 52°. Head motion was minimized with foam padding.

During MRI scanning, subjects were supine inside the magnet. Visual stimuli were presented through an MRI compatible eyepiece to the right eye of the subject. Stimuli were delivered from the digital light synthesizer through a randomized optical fiber cable (Fiberoptics Technology Inc). The randomization of the fiber optic cable helped to minimize any spatial inhomogeneities in the spectrum of the stimulus. The eye piece provided adjustable focus to account for variation in refractive error. As the stimulus was a spatially uniform field, however, the effect of any spatial blur upon the stimulus was minimal.

Subjects used either button of a two button MR compatible response device (Current Designs) to respond to attention events during the functional runs.

Both anatomical and functional data were preprocessed according to the HCP minimal preprocessing pipeline (Glasser et al., 2013) Briefly, the anatomical data were passed through the pre-freesurfer, freesurfer, and post-freesurfer steps of the HCP minimal preprocessing pipeline. This was used to create an MNI registration, a Freesurfer segmentation (Dale et al., 1999; Fischl et al., 1999), and a surface mesh. The functional data were preprocessed with both the volume and surface pipelines. The volume pipeline applied gradient distortion correction, motion correction (FLIRT 6 DoF; Jenkinson et al., 2002), top-up phase encoding distortion correction (Smith et al., 2004), and registered the functional images to the anatomical images. The surface pipeline mapped the volume data to the CIFTI grayordinate space which includes a left and right 32K cortical surface mesh and subcortical voxels. Finally, the functional data were passed through the ICAFIX pipeline, which uses independent component analysis and noise/not-noise classification to denoise the time course.

After initial preprocessing, we performed a series of subsequent steps before analyzing the time course data. We used a V1 region of interest (ROI) to extract the time series from primary visual cortex (see below for definition of retinotopic maps). The signals from each voxel were mean centered and converted to percent signal change. We then performed nuisance regression using the relative motion estimates and attentional events as regressors. The relative motion regressors (translation of X, Y, and Z and yaw, pitch, and roll) were mean centered and scaled by their respective standard deviations. The attention event regressor was modeled as a series of delta functions located within the TRs in which the events occurred, convolved with a hemodynamic response function. The nuisance regression was performed using the MATLAB linear regression function mldivide (MathWorks) with the residual of the model used as the ‘cleaned’ timed series.

Next, we time-point censored the time series of all voxels based on the motion estimates. This was done using a modified version of Power et al., 2014. We converted the relative yaw, pitch, and roll estimates to millimeters by using the distance that each of these angles subtend on an assumed 50 mm sphere. We then took the L2-norm of the six translation and mm rotation estimates as a metric of frame-wise displacement (FD). We censored three contiguous TRs centered on any time point with an FD > 0.5. Time points that exceeded the threshold were excluded from analysis, and a table of the number of censored frames can be found in Table 4. Finally, we applied polynomial detrending by fitting a fifth order polynomial to the time course from each voxel and subtracting it from the signal. Analyses performed at the level of V1 were done using the median value across voxels at each time point to represent the V1 signal. Vertex-wise analyses were performed on the preprocessed time course of individual vertices.

Retinotopic regions of interest (ROI) were defined using the anatomical template neuropythy (Benson et al., 2014) which provides eccentricity and polar angle maps for V1, V2, and V3. From this atlas, we defined a V1 ROI using the voxels in area V1 between 0–20° eccentricity and 0–180° polar angle. We set the eccentricity upper bound of the ROI to be 20° to provide a conservative boundary to ensure that we only analyzed stimulated vertices. This accounts for some uncertainty in the exact retinal size of the stimulated area due to, for example, variation in the distance of the eyepiece to the eye of the subject.

To define the retinotopic regions of interest used for hV4/Vo1 analysis, we registered the Wang et al., 2015 retinotopic atlas to CIFTI space. This was achieved using the subject specific atlas as defined by Neuropythy (Benson and Winawer, 2018). This atlas was then registered to MNI space through the use of ANTs (Avants et al., 2008). With the atlas in MNI space, we then transformed it to CITFI greyordinate space and used the hV4 and the VO1 ROIs to extract the time series data.

We derived subject-specific hemodynamic response functions (HRFs) for each subject and measurement set. The HRF was derived using all the functional runs within a measurement set, using the V1 region of interest median data. The time-series data were fit with a non-linear model that simultaneously estimated the beta weights of the GLM to fit stimulus responses, and the parameters of an HRF model. The HRF model was composed of the first three components of the ‘FLOBS’ basis set (Woolrich et al., 2004). The best fitting HRF model was then used to fit that subject’s data for both the GLM and QCM models.

We used an ordinary least squares regression with a stimulus design matrix that described the stimulus order of a run. The regression matrix contained one regressor per stimulus block as well as a single regressor for the baseline, with the length of the regressor equal to the number of timepoints (360 TRs). The regressor for each stimulus block in a run was modeled by a binary vector that indicated the timepoints when the stimulus was present (ones) or absent (zeros), convolved with the HRF. This resulted in 21 GLM beta weights per run. For all analyses, model fitting was performed using the concatenation of all functional runs within a measurement set. For fitting of the GLM, we concatenated the stimulus design matrices such that like contrast/direction pairs were modeled by the same regressor. The concatenated stimulus design matrix for a measurement set had a total of 41 regressors (20 direction/contrast pairs from Session 1, 20 direction/contrast pairs from Session 2, and a shared baseline regressor) and 7200 timepoints. Additional details of the GLM are provided in the Appendix 1 section.

To obtain contrast response function for each color direction, the time course data for each run was fit using a general linear model to obtain the effect that each stimulus had on the BOLD fMRI response. Grouping the GLM beta weights by chromatic direction defined a set of 8 contrast response functions (CRFs), one per direction. A CRF describes the relationship, within a particular chromatic direction, between the contrast of the stimulus and the measured response. The CRFs obtained using the GLM beta weights fit to the concatenated time series of Subject 2 can be seen in Figure 3. The panels of Figure 3 show the CRFs for the eight different chromatic directions. The CRFs for Subject 1 and 2 can be seen in Figure 2—figure supplements 1–2.

The error bars and error regions in all figures are the 68 percent confidence intervals. We used 68% confidence intervals as this approximates +/- 1 SEM for a normal distribution. The percentiles we use to estimate error are the result of a bootstrap analysis. The bootstrap analysis was implemented as random sampling with replacement of the runs within a measurement set. The randomly drawn runs in both sessions were concatenated and fit by all models. We performed 200 bootstrap iterations and identified the 68% percent confidence interval from the bootstrap results.

To evaluate model performance, we employed a leave-runs-out cross-validation strategy. For each cross-validation iteration, runs from Session 1 and Session 2 were randomly paired within the same measurement set. These pairs of runs were held out and the models were fit to the remaining 18 runs. From these model fits, a time course prediction for the left-out runs were obtained from both models. We computed the R^2 value between these predictions and the time course of the held-out runs. The average R^2 value across the 10 cross-validation iterations was used to compare models.

To evaluate the generalizability of the QCM, we implemented a leave-session-out cross-validation. As the eight chromatic directions tested were separated into two sessions with the same grouping across all subjects and measurement sets, we could evaluate the ability of the QCM to make predictions for chromatic directions and contrasts not used to fit the model parameters. Within each measurement set, we fit the QCM to Sessions 1 and 2 separately and evaluated how well the parameters of the model predicted responses to stimulus directions in the held-out session. These predicted responses were grouped by chromatic direction in order to construct a set of contrast response functions. The error bars in the CRFs are the 68% confidence intervals computed using bootstrapping, where we randomly sampled runs with replacement and compute the leave-session-out cross-validation a total a 200 times.

Cortical surface maps were generated to visualize the ellipse angle and minor axis ratio parameters of the QCM on the V1 cortical surface. To generate the surface parameter maps, we fit the QCM, within a single vertex, to the concatenated time course from all runs in a measurement set. We repeat this fit for all vertices within the visual areas template map from neuropythy (Benson et al., 2014). To visualize the surface parameter maps in a manner that highlights differences in fits as a function of cortical position, we created a series of scatter plots that relate the minor axis ratio and angle parameters to the eccentricity of their respective vertices. The regression lines in the scatter plots are robust regression lines implemented through the built in MATLAB function robustfit which adaptively reweights the data to discount the effects of outliers.

We fit QCM the model to the concatenated time series for both sessions within a measurement set for each subject. The data were fit using the MATLAB function fmincon to find a set of model parameters that minimize the difference between the actual fMRI time course and the QCM prediction of the time course, computed as the root mean squared error.

All code used to perform analyses in this paper may be found in our public GitHub repository: https://github.com/BrainardLab/LFContrastAnalysis_eLife.git (copy archived at swh:1:rev:c9b4cbe72e69e4d3d623ec8b7fc62076e2fe1a22, Barnett, 2021).

Tables 5–7 show the stimulus validation measurements for all subjects and sessions. The tables show the mean and standard deviation of stimulus vector angles and lengths computed from 10 validation measurements (five pre-experiment, five post-experiment). Center and periphery denote which set of cone fundamentals were used to calculate cone contrast of the stimuli referring either the 2° or 15° CIE fundamentals, respectively.

The model used to provide a benchmark for the quadratic color model (QCM) is the general linear model (GLM). The GLM has the form:

This states that the measurement ($Y$) is equal to the model matrix ($X$) times the weights ($\beta$) plus the residual error ($ϵ$). For our GLM, $Y$ is a column vector of the concatenated time series from all the fMRI runs within a measurement set (20 runs):

The subscript $i$ indicates a particular time point, with the subscript $t$ being the total number of time points in a measurement set. Here a time point corresponds to one TR of the BOLD response, and t = 7200 (20 runs with 360 TRs). Our model comparison was for the aggregated V1 response, and we took each element y_i of $Y$ to be the median BOLD fMRI response for the corresponding TR, with the median taken across the voxels in the V1 ROI.

The model matrix $X$ contains the predictor variables for the linear model. Each column of $X$ is a regressor corresponding to a single stimulus (a single chromatic direction/contrast pair or the baseline (0 contrast) uniform field). The regressors are created by convolving a binary indicator vector with the hemodynamic response function (HRF) that accounts for the sluggish BOLD response to a stimulus event. The binary indicator vectors contain 1 when the stimulus was present and 0 otherwise. The convolution then produces a predictor of the measured BOLD fMRI response for that stimulus condition. In our study, $X$ has 41 columns corresponding to the pairing of the eight chromatic directions and five contrasts levels plus the baseline

The general linear model predicts $Y$ as a weighted linear combination of the regressors in the columns of $X$. The weight applied to each regressor is given by the elements of $\beta$:

Using $\beta$ and $X$, we can generate a predicted time course $\widehat{Y}$. We take each ${\beta }_i$ to be a proxy for the aggregate V1 BOLD fMRI response for the corresponding stimulus condition. This interpretation is associated with our particular choice of scaling the indicator variables in the regressors before convolution with the HRF (that is the choice of 1 for the TRs during which the stimulus was present). We determined $\beta$ using linear regression as implemented in the MATLAB mldivide operation, $\beta =X\backslash Y$ (see https://www.mathworks.com/help/matlab/ref/mldivide.html). The GLM prediction of the BOLD fMRI response at each time point y_i is therefore:

The quadratic color model allows for the prediction of the BOLD fMRI response to any stimulus that lies in the LM contrast plane. This includes predictions for stimuli not in the measurement set used to fit the model. The QCM makes its predictions through three steps. The first step calculates the ‘equivalent contrast’ of the stimulus which can be thought of as the effective contrast of the stimulus in V1 accounting for the differences in sensitivity across chromatic directions. Equivalent contrast is the vector length of the stimuli after this linear transformation. The next step is a response non-linearity applied to the equivalent contrast to predict the underlying neural response. Finally, a convolution of this underlying response with the HRF results in predictions of the BOLD fMRI response. Here, we provide explanations and equations for the model.

We start by considering a stimulus modulation whose predicted BOLD fMRI response we wish to know. A stimulus modulation is denoted by the column vector $c$ whose two entries are the L and the M cone contrast of the stimulus ($l$ and $m$, respectively): 

The figure below outlines how we transform such a stimulus to a response, prior to convolution with the HRF. Panel A shows a 3-dimensional contrast-response space with the ($x$,$y$) plane representing the LM contrast plane and the $z$ axis giving the response $r$ corresponding to each point in the LM contrast plane. In this representation, all possible stimuli and responses form an inverted bell shape surface, as illustrated below. At constant values of $r$ (constant height on the $z$ axis), cross sections through this surface shows the elliptical isoresponse contours of the QCM. Each elliptical isoresponse contour describes the set of LM contrast combinations that elicit the same response $r$. The teal and red dots shown in the LM plane represent example stimulus modulations, chosen in two color directions at contrasts corresponding to the five (dark blue, blue, aqua, green, yellow) isoresponse contours illustrated. 

To obtain the equivalent contrast corresponding to stimulus $c$, we first apply a linear transformation $M$ to $c$ to obtain a transformed representation of the stimulus, $e=Mc$. The vector $e$ is a two-dimensional column vector whose entries we refer to as e₁ and e₂. We call the transformed representation of the stimulus the equivalent contrast space, and as we show below we choose the linear transformation $M$ such that in this space the isoresponse contours are circles. Panel B of the figure shows the same information as in Panel A represented in the equivalent contrast space. Here the ($x$,$y$) plane gives the values of e₁ and e₂, and the isoresponse contours plotted with respect to the equivalent contrast plane are circular. Note that in the equivalent contrast plane, the distances between the plotted teal points are the same as the distances between the corresponding plotted red points. This is a consequence of the fact that the transformation $M$ is chosen to make the isoresponse contours circular.

Appendix 1—figure 1

Download asset Open asset

We define the equivalent contrast $k$ of a stimulus as the vector length of the transformed stimulus $e$, $k=|e|$. 

The next step of the model is a non-linearity that maps between equivalent contrast ($k$) and BOLD fMRI response (dashed black line in Panel C of the figure). This is possible since all stimuli in equivalent contrast space with the same equivalent contrast predict the same underlying response, regardless of the chromatic direction of the stimuli. Therefore, we can focus solely on the relationship between equivalent contrast and the associated response. We call the static non-linearity the ‘equivalent contrast-response function’. The teal and red closed circles shown in Panel C correspond to both the teal and red points shown in Panels A and B, and these overlap since they lie on the same set of isoresponse contours. 

The linear transformation $M$ needed to compute the equivalent contrast representation of $c$ is derived by starting with the equation for an ellipse centered at the origin, in which a positive-definite quadratic function of the cone contrasts is equal to a constant for all points on the ellipse. This equation is given below (left hand side) and re-expressed in matrix-vector form (right hand side).

Changing the value of the constant $k^2$ changes the scale of the ellipse without changing its shape, and as we will see below $k$ is the equivalent contrast corresponding to each constant-shape elliptical isoresponse contour. We rewrite the matrix representation of the elliptical locus as: 

This yields

Because the coefficients $A$, $B$, and $C$ are constrained so that $Q$ is a symmetric positive definite matrix, $Q$ can be factored through its eigenvalue decomposition and rewritten as

where $V$ is an orthonormal (rotation) matrix and $\mathrm{\Lambda }$ is a diagonal matrix with positive entries. Since $\mathrm{\Lambda }$ is a diagonal matrix we can further decompose as

where $S$ is diagonal with entries equal to the square root of the corresponding entries of $\mathrm{\Lambda }$. 

The matrix $V$ expresses a rotation in the cone contrast plane and may be parameterized by the rotation angle p₁:

The matrix $S$ is diagonal and when applied to a vector simply scales the entries of that vector. Although $S$ normally has two degrees of freedom in the general case, we lock the scale of the major axis to 1. Therefore we are only concerned with the scaling of the minor axis (p₂). Thus we can define

In this formulation, p₁ and p₂ parameterize the shape of an elliptical isoresponse contour and $k$ parameterizes its scale. The parameter p₁ is what we call the angle of the major axis in the main text, while p₂ is the minor axis ratio. 

We set $M=S^T{V}^T$ and rewrite $Q$ as:

The matrix $M$ then transforms the cone contrast vector $c$ to the equivalent contrast vector through $e=Mc$. Recall that equivalent contrast is the vector length of $e$. Therefore, the equivalent contrast of the points on the ellipse corresponding to $k$ is given by

To see this, note that

Thus, given $Q$, we can compute the equivalent contrast for any stimulus $c$ and apply the equivalent contrast-response function to predict the its response. The specific non-linearity we use for the contrast-response function is a Naka-Rushton function. This function is a four parameter saturating non-linearity that is defined by the following: 

The neural response ($r$) is a function of the equivalent contrast ($k$). The parameters of the Naka-Rushton function are the amplitude ($a$), exponent ($n$), semi-saturation ($s$), and the offset ($h$). These parameters control the gain, the slope, the position along the x axis, and the y-axis offset of the non-linearity, respectively. The shape of the non-linearity controls how the underlying response changes with equivalent contrast. 

Finally, we need to convert the neural response to a prediction of the BOLD fMRI response($\widehat{Y}$). To do this, we convolve the underlying response with the hemodynamic response function (HRF).

Predictions of the BOLD fMRI response via the QCM are thus made using six parameters: angle (p₁), minor axis ratio (p₂), amplitude ($p_3=a$), exponent ($p_4=n$), semi-saturation ($p_5=s$), and offset ($p_6=h$). We can define a parameter vector $P$ for the QCM as:

We fit $P$ to the data by using a non-linear parameter search routine fmincon (Matlab, see https://www.mathworks.com/help/optim/ug/fmincon.html) to find the parameter vector that that minimizes the difference between the measured BOLD fMRI time course $Y$ and the prediction $\widehat{Y}$ obtained using the QCM. This takes the form of:

where the objective function we are minimizing is the root mean squared error. Once the values of ${\displaystyle P^{\ast }}$ are found for an individual subject, we can use them to predict the BOLD fMRI response to any $c$ within the LM cone contrast plane.

The linear channels model (LCM) is based on the work of Brouwer and Heeger, 2009. They develop a model that decodes the hue angle of a stimulus a subject was viewing from BOLD fMRI measurements. The stimuli in their study were modulations in the isoluminant plane of the CIELAB color space, and were parameterized by hue angle within this plane. Their model is based on linear channels tuned to the this angular representation, and it did not explicitly consider the effect of stimulus contrast. This was sufficient for their analysis, as they studied only a single contrast for each angle. To develop a version of the LCM that may be applied to our stimuli, we need both to adopt the concepts to apply to stimuli in the LM contrast plane and to add a model component that handles contrast.

The angles used in our implementation of the LCM are angles in the LM contrast plane, with 0 corresponding to the positive abscissa, rather than angles in the CIELAB isoluminant plane. These angles are computed from the L- and M-cone contrasts of our stimulus modulations, as the arctangent of the ratio of M- to L-cone contrast:

where $\theta$ is then the angle corresponding to one of our stimulus modulations. 

Given the angular stimulus representation, each linear channel is characterized by its sensitivity to stimulation from each possible stimulus angle - basically a tuning function over angle. More specifically, the tuning functions are implemented as half-rectified cosines raised to an exponent. What differs across each channel is the phase of the tuning function (the peak sensitivity). Thus the sensitivity of the $i^{th}$ channel may be written as: 

where $i$ indicates the channel, $n$ is the exponent which controls the narrowness of the channel tuning functions, and ${\varphi }_i$ is the angle of peak sensitivity of the $i^{th}$ channel. In calculations, the angle $\theta$ is discretized, and both stimuli and mechanism tuning can be represented in matrix-vector form, as we describe below.

The first step of the model is to compute the channel responses. This is done with a vector representation of the stimulus. Each entry of the stimulus row vector represents one angle, and for the stimulus modulations we used, only one entry is non-zero for any given modulation. In our computations, we discretized angle in one degree steps (the resolution of the representation of $\theta$). The magnitude of the non-zero entry represents the contrast of the modulation. Similarly, the sensitivity of a channel may be represented by a column vector over the same angular discretization, with each entry of the vector being the sensitivity of the channel at the corresponding angle. The dot product of a stimulus vector with a channel vector yields the response of the channel to the stimulus.

The set of channels can be represented by the columns of a matrix C (nAngles x nChannels) and the set of stimuli represented by a matrix S (nStimuli x nAngles). Therefore, the hypothetical channel outputs, $H$ (nStimuli x nChannels), can be calculated as:

The overall LCM response to a given stimulus is given as a weighted sum of the individual channel responses ($H$), with the weights, $w$ (nChannels x nVoxels) acting as parameters of the model. If we consider responses across a set of voxels, as was done by Brouwer and Heeger, 2009, the response $b$ (nStimuli x nVoxels) for a set of stimuli is given by: 

In our analysis, we fit the LCM to the median time course of V1 and therefore we set $nVoxels=1$.

In the context of the LCM, isoresponse contours are made up by the set of stimuli that satisfy the following:

where $k$ is a constant target response. 

In their work, Brouwer and Heeger, 2009 used six channels and an exponent of 2. The channels had peak sensitivities at 0°, 60°, 120°, 180°, 240° and 300°. Since our stimulus modulations were symmetric around around a background, we enforce that channels offset by 180° used the same weights in the linear combination. This results in three symmetric channels with pairs located at 0° and 180°, 60° and 240°, 120° and 300°.

One key difference between our experiment and the the original Brouwer and Heeger, 2009 work is that we varied the stimulus contrast in each chromatic direction. To extend the LCM to handle contrast, we treated the overall LCM response as an equivalent contrast and passed this through a common Naka-Rushton function. This approach mirrors how we handled contrast in the the QCM. To model the BOLD fMRI response, we convolve the output of the Naka-Rushton with the HRF. We fit the LCM to our data using a method analogous to the one used to fit the QCM. More specifically, we found the three linear channel weights and the Naka-Rushton function parameters that yielded the best fit to the BOLD response time course. These parameters were found simultaneously through the use of MATLAB’s fmincon optimization routine. We also fit a variant of the LCM analogous to a model used by Kim et al., 2020 with sharper ($co{s}^6$) channel tuning and using eight rather than six underlying mechanisms. 

Both versions of the LCM yield isoresponse contours similar in shape to those obtained with the QCM (Figure 6—figure supplement 1), with cross-validated $R^2$ values essentially the same as those obtained using the QCM (Figure 4—figure supplement 1). Notably, that the version of the LCM with more channels yielded an isoresponse contour that more closely approximated the isoresponse contour of the QCM than did the version of the LCM with fewer channels. 

A key difference between the LCM and the QCM is in how the functional properties of the model relate to the L- and M-cone contrasts that characterize the early visual system responses to the stimuli. As discussed in the main text, the QCM can be implemented as the sum of the squared responses of two mechanisms, each of which responds as weighted sum of L- and M-cone contrast. Thus underlying the QCM is a quadratic (squaring) non-linearity. In the LCM, on the other hand, the underlying linear channels are tuned for angle in the LM contrast plane. Since angle is obtained as the arctangent of the ratio of M- and L-cone contrasts, the LCM is based on a different form of non-linearity than the QCM, which would entail different neural computations. To put it another way, the apparently simple linear form of the LCM when expressed in terms of stimulus angle is less simple when referred back to the L- and M-cone contrast representation. This same basic point applies to LCMs formulated with respect to hue angle in the CIELAB isoluminant plane, as in Brouwer and Heeger, 2009 and Kim et al., 2020. Measurements of the response to mixtures of modulations might be used in the future to differentiate between the non-linearities embodied by the LCM and the QCM.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This work has been supported by the National Institutes of Health (Grants:RO1 EY10016 and Core GrantP30 EY001583; https://www.nih.gov/) and National Science Foundation Graduate Research Fellowship (DGE-1845298). We thank Joris Vincent, Jack Ryan, Ozenc Taskin, and Nicolas Cottaris for technical assistance.

Human subjects: The research was approved by the University of Pennsylvania Institutional Review Board (Protocol: Photoreceptor directed light modulation 817774). All subjects gave informed written consent and were financially compensated for their participation.

1. Preprint posted: December 4, 2020 (view preprint)
2. Received: December 9, 2020
3. Accepted: July 27, 2021
4. Accepted Manuscript published: August 3, 2021 (version 1)
5. Version of Record published: September 20, 2021 (version 2)

© 2021, Barnett et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.