Modeling neural contrast sensitivity functions in human visual cortex

The detection of contrast and spatial frequency are fundamental aspects of vision. These features are important for object recognition and essential in everyday life. Where contrast is important in distinguishing an object from its background, spatial frequency reflects object size (Chung & Legge, 2016). The ability to detect a certain contrast is given by contrast sensitivity (100/contrast threshold) and depends on spatial frequency (Campbell & Robson, 1968; Lesmes et al., 2010; Sowden et al., 2002). The relation between contrast sensitivity and spatial frequency is described by the contrast sensitivity function (CSF). Behaviorally, the CSF defines the lowest contrast one can perceive as a function of spatial frequency (Campbell & Robson, 1968), thereby defining the threshold between the visible and invisible (Pelli & Bex, 2013).

The CSF is used to assess visual function. For example, the CSF is affected in many ophthalmological conditions, including amblyopia (Howell et al., 1983; Koskela & Hyvarinen, 1986; Sjöstrand, 1981; Wang et al., 2017), macular degeneration (Kleiner et al., 1988), optic neuritis (Zimmern et al., 1979), glaucoma (Ichhpujani et al., 2020), retinitis pigmentosa (Hyvärinen, 1983), cataract (Vasavada et al., 2014), and corneal edema (Hess & Garner, 1977). Additionally, the CSF can be altered in neurological conditions such as multiple sclerosis (Regan et al., 1981), cerebral lesions (Milling et al., 2014), Parkinson’s disease (Ridder et al., 2017), and schizophrenia (Cimmer et al., 2006).

In general, the CSF correlates with visual acuity (Hou et al., 2010; Stalin & Dalton, 2020). However, the CSF measures the detection of a stimulus using a wide range of contrasts and spatial frequencies, and has therefore been proposed as a more suitable tool to assess visual performance compared to standard visual acuity tests (Huang et al., 2007). For example, the CSF measures visual deficits missed by standard visual acuity tests. Thus, the CSF provides deeper insights regarding functional vision and improves detection of visual pathology (Huang et al., 2007; Lesmes et al., 2010).

The CSF characterizes visual perception and is altered by changes in the eye as well as changes in neural processing. In healthy participants, the CSF can be altered by cognitive manipulations such as attention, for example, contrast sensitivity is increased for attended stimuli and decreased for unattended stimuli (Pestilli et al., 2007). Furthermore, the CSF changes across development, during infancy as well as late childhood (Dekker et al., 2020). In addition, several visual disorders are at least, in part, caused by neural deficits, for example amblyopia (Barrett et al., 2004) and glaucoma (Murphy et al., 2016). This highlights the relevance of studying neural processes underlying the CSF, in addition to the ocular components. Thus, a neural measure of the CSF would expand our understanding of how cortical processing is altered with cognition and different visual disorders.

Here, we translate the CSF from psychophysics, by introducing a new method that estimates the neural CSF (nCSF) in the human visual cortex using fMRI. The nCSF concept is derived from ophthalmology where it refers to the ability of the retina together with the brain to resolve an image, that is, not affected by eye-optics (Campbell & Gubisch, 1966; Le Grand, 1935; Michael et al., 2011). The nCSF approach builds on previous studies showing that individual neurons in the visual cortex are sensitive to contrast and spatial frequency (Albrecht & Hamilton, 1982; Levitt et al., 1994; Sclar et al., 1990). This sensitivity is also present at the neural population level with fMRI studies showing systematic changes in contrast (Boynton et al., 1999; Marquardt et al., 2018), spatial frequency (Aghajari et al., 2020; Broderick et al., 2022; Henriksson et al., 2008; Singh et al., 2000; Sirovich & Uglesich, 2004), and their combination (Goulet & Farivar, 2024). Here, we combine population measures of the CSF (akin to the population receptive field (pRF) method; Dumoulin & Wandell, 2008) and we model the transition from no response to full response, by varying the slope of the contrast response function (CRF, Boynton et al., 1999).

We show that our nCSF model effectively captures the variance in the fMRI time series. Moreover, the properties of the nCSF vary systematically with eccentricity, polar angle and across the visual hierarchy. Overall, we describe a quantitative validated framework to model nCSF properties.

We present data from five participants (two males, age range 26-45 years). All participants had normal or corrected-to-normal visual acuity, as confirmed by testing visual acuity using a tumbling E eye chart. The participants gave written informed consent prior to the start of the experiment. The study was approved by the Ethical Committee of Vrije Universiteit Amsterdam in accordance with the World Medical Association’s Declaration of Helsinki.

The stimuli were generated using PsychToolbox (Brainard, 1997; Pelli, 1997) in MATLAB (version R2018b, Mathworks). Participants viewed the stimuli on a gamma corrected 32-inch BOLD screen display (Cambridge Research, 1920 x 1080 pixels), through an angled mirror attached to the head coil (distance 220 cm). 10-bit mode was enabled on the BOLD screen display to facilitate presentation of stimuli with contrasts below 0.78% Michelson contrast. Participants were instructed to fixate on a dot in the center of the screen, and to press a button when the dot changed color.

The nCSF stimulus consisted of different static sinewave gratings with a circular aperture of 10.2 degree diameter of visual angle and raised cosine edge, presented on a mean luminance background. The gratings varied systematically in contrast and spatial frequency; see Figure 1. Each stimulus block lasted 18 s and consisted of gratings with the same spatial frequency but systematically changing (either decreasing or increasing) in contrast (Fig. 1A). Three different gratings were shown every 1.5 s with the same contrast but a different orientation and phase (Fig. 1B), to ensure that the stimulus was refreshed on the retina. Twelve different contrasts were used in each stimulus block, ranging from 0.25–80% Michelson contrast. An increasing and decreasing contrast block was included for each spatial frequency (0.5, 1, 3, 6, 12, 18 c/deg), and the order in which these blocks were presented was fixed; see Figure 1A. We included stimulus blocks with both increasing and decreasing levels of contrast to avoid directional effects introduced by the HRF (similar as in pRF designs), and to avoid the threshold being dependent on whether the contrast change is increasing or decreasing.

Before every two stimulus blocks there was a mean luminance block lasting for 15 s, and there was an extra 15 s mean luminance block in the end. The orientation of the gratings was randomized, and two subsequent gratings had at least a 45-degree difference in orientation. Each grating presentation lasted for 300 ms and was followed by a 200 ms presentation of mean luminance (see Fig. 1B). The total stimulus protocol lasted 321 s per run. The contrast range presented was dependent on the spatial frequency shown in the stimulus block, ensuring an optimal sampling of the CSF (see Fig. 1C, and Table S1).

Retinotopic (pRF) data were collected as part of previous experiments; the stimuli consisted of a drifting bar, moving in eight directions. For four of the participants, the texture inside the bar comprised natural images. For the fifth participant, the bar contained flashing checkerboards. Stimuli were presented using the same screen (without the 10-bit mode) and software as in the nCSF experiment.

Anatomical and functional magnetic resonance imaging (fMRI) data were acquired on a 7 Tesla Philips Achieva scanner (Philips, Best, Netherlands) with an 8-channel MultiX head coil. Anatomical images were collected using a 3D MP2RAGE sequence (TR = 6.2 ms, TE = 3.0 ms, flip angle = 5.0 degrees, FOV = 220 x 220 x 164 mm, voxel size = 0.7 mm isotropic) (Marques et al., 2010).

fMRI data were acquired using a T2*-weighted 2D echo-planar imaging (EPI) sequence (TR = 1.5 s, TE = 22.5 ms, flip angle = 65 degrees, FOV = 216 x 216 x 97 mm, 57 slices, voxel size = 1.7 mm isotropic) oriented across the visual cortex. For each participant, 10–12 functional scans were acquired, with an approximate duration of 5 min per functional run. Each functional run was followed by a TOPUP scan in order to correct for local image distortions. Additionally, pRF data were collected for one participant (8 functional runs) with the same sequence for visual field mapping (see Section 2.6). For the remaining four participants, pRF mapping data were already available.

For the anatomical data, the MP2RAGE sequence was used to obtain a T1-weighted anatomical image, resulting in a set of different gradient echo images, that is, the first inversion (T1-weighted image) and a second inversion (proton density scan). The first and second inversion were combined to generate a single anatomical image, corrected for proton density. The anatomical image was skull-stripped using the 3dSkullstrip function in AFNI (Cox, 1996). Segmentation of gray and white matter was performed using CBS Tools in MIPAV (Bazin et al., 2007). The cortical surfaces for each participant were reconstructed using the FreeSurfer 7.2 recon-all function (Dale, 1999). To improve the FreeSurfer reconstructions particularly around the sinus, segmentations from CBS Tools were added to Freesurfer’s “brainmask.mgz” file, and recon-all was run again.

For the functional data, all functional images were pre-processed using AFNI. The warp field of all functional images was estimated based on the EPI and TOPUP of each run to correct for susceptibility distortions. The motion parameters were calculated to apply motion correction to all functional images. Then, all warped and motion corrected volumes were combined to generate an average EPI image. This EPI image was registered to the anatomical image. The volume of the mean EPI image was masked using the AFNI function 3dAutomask to reduce volume size, and was zero-padded using 3dZeroPad. The center of mass of the anatomical volume was aligned to the mean EPI image using @Align_centers to improve co-registration. Then, the mean EPI image was manually shifted and rotated using the Nudge dataset AFNI plugin, after which the automated registration function 3drotate optimized the co-registration using affine transformation. The co-registration was then applied to all EPI images of all individual functional runs using the function 3dNwarpApply with nearest-neighbor interpolation. The co-registered functional volumes were then projected to the participant’s cortical surface using the mri_vol2surf function in FreeSurfer (Dale, 1999). Veins were identified by displaying the mean EPI signal on the cortical surface, and masking those vertices with relatively low signal intensity. The BOLD time series were then detrended by demeaning, applying the discrete cosine transform (DCT), and removing the first three DCT coefficients to eliminate low-frequency trends. The detrended data were further processed by averaging over functional runs and converting to BOLD percent signal change.

The nCSF model combines the contrast sensitivity function (CSF) with the contrast response function (CRF) and aims to fit the CSF in each cortical location (voxel), thereby obtaining the nCSF (Fig. 2). The CSF is described by Chung and Legge (2016) and is characterized as an asymmetric parabolic function with the following equation:

where f(SF) is the contrast sensitivity at spatial frequency (SF), CS_p is the peak contrast sensitivity, SF_p is the spatial frequency at which CS_p occurs (peak spatial frequency), and width_L and width_R are the curvatures of the left and right branches of the asymmetric parabolic function, respectively (Chung & Legge, 2016). Note, following the notation in Chung and Legge (2016), “width” parameters refer to the “fall off rate”, hence a CSF with a high width_r actually has a narrower shape (there is a faster drop in sensitivity to higher spatial frequencies). The width_L is set at 0.68 (Chung & Legge, 2016), and width_R is varied as a model parameter, since we expect the most variance on the right branch at higher spatial frequencies (SF > SF_p). We report the normalized area under the log contrast sensitivity function (AUC, %) as a summary metric of the full CSF for spatial frequencies between 0.5–18 c/deg (stimulus space). AUC is calculated by approximating the integral of the logarithmic CSF, and normalizing it with reference to the area under a “standard” logarithmic CSF calculated from healthy controls (from Chung & Legge, 2016: CS_p = 166, SF_p = 2.5 (c/deg), width_L = 0.68, width_R = 1.28, shown in Fig. 1C).

The parameters SF_p and CS_p are transformed to log₁₀ space before fitting the CSF (Equation 1, Chung & Legge, 2016). Here, we used an adapted version of Equation 1 where all parameters are defined in linear space before fitting the CSF:

For transforming the binary response of the CSF to a gradual response, we added the CRF to the nCSF model (Fig. 2). We fitted CRFs using the following equation (modified from Boynton et al. (1999)):

where R is the fMRI response and C is the amount of RMS-contrast. The variables for Q and q define the shape of the CRF. Q represents the contrast where the fMRI response is at 50%, and was depended on the CSF. We fitted the variable q or slope_CRF.

The nCSF parameters were estimated from the fMRI data using a model-based fitting approach (Fig. 2), similar to the population receptive field (pRF) method (Dumoulin & Wandell, 2008). First, the fMRI blood oxygen level dependent (BOLD) response is predicted through multiplying the nCSF model with the stimulus sequence and convolving the time course with the hemodynamic response function (HRF). Here, the HRF was modeled using the two gamma basis functions (Glover, 1999). The coefficients corresponding to the canonical HRF, its derivative, and dispersion were fixed at values of 1, 1, and 0, respectively (Pedregosa et al., 2015). The same HRF was used in fitting both the nCSF and pRF model. We varied the model parameters (CS_p, SF_p, width_R and slope_CRF) in a coarse-to-fine manner. All model parameters were varied until the residual sum-of-squares (RSS) is minimized between the predicted and measured fMRI data, thereby maximizing the variance explained (r²) of the model. The nCSF model was fitted to all cortical locations (voxels) independently. We selected cortical locations with r² > 30% in both pRF and nCSF model fits for further analyses.

For the region of interest (ROI) definition, pRF mapping was used (Dumoulin & Wandell, 2008). These data were already available for four participants, and we collected pRF data for one participant (see Section 2.3). We included visual field maps V1, V2, V3, hV4 and visual field map clusters TO, LO, and V3AB as ROIs (Wandell et al., 2007).

We validated the method using simulations (Dumoulin & Wandell, 2008; Lerma-Usabiaga et al., 2020), where the ground truth is known, and estimated the accuracy through model parameter recovery. To validate the nCSF method, we simulated two nCSF models with different CSF parameters (in green: SF_p = 1 c/deg, CS_p = 150, width_R = 1.3, in red: SF_p = 2 c/deg, CS_p = 100, width_R = 1). For each simulated dataset, the predicted response was calculated using the stimulus sequence (as described in Section 2.2) and combination of model parameters, convolved with the standard HRF. We normalized these synthetic voxel time series and added three different noise levels. Noise was generated by sampling from a normal distribution (M = 0, SD = 1) multiplied by a scaling factor 0.4, 0.7, or 1.1. For each noise level and combination of model parameters, 100 time series were generated. Next, we fitted the nCSF model using the same HRF used for the main data analysis (see Section 2.5). The model parameters were extracted from the simulated datasets, and the variance explained for each of the three noise levels was computed, resulting in the low, medium, and high variance explained categories. The two different nCSF curves for the three variance explained categories are shown in Figure 3.

Additionally, we simulated the effect of the HRF on the model parameters (see HRF Simulation section in Supplementary Materials). Furthermore, we simulated how slope_CRF varies with noise, and the effect of slowing speed at which the stimulus changes contrast (see CRF Simulation section in Supplementary Materials).

Validation is essential for method development (Aqil et al., 2021; Dumoulin et al., 2003; Dumoulin & Wandell, 2008). Modern data-analyses techniques are complex and it is impossible to check the software by eye: ground-truth datasets are needed (Lerma-Usabiaga et al., 2020). Here, we validated the nCSF model using simulations.

We created two “ground-truth” nCSF models, with different parameters (red and green in Fig. 3), added noise, and refitted the resulting time series to determine how well we can recapture the model parameters (see Section 2). We splitted the resulting fits into three variance explained levels: high (r² > 50%), medium (30% < r² < 50%), and low (10% < r² < 30%). Across all variance explained levels (see Fig. 3A-C), the median estimated nCSF curves (dashed lines in red and green) closely match the true curve (solid lines in red and green). However, the spread of the estimated nCSF curves increased with lower levels of variance explained (compare shaded regions in Fig. 3A-C). The median values of the recovered parameters closely match the ground-truth values (compare the dotted with solid horizontal lines in Fig. 3D-I). Again, the spread of the recovered parameters increases from high to low variance explained levels (see Fig. 3D-I).

We fitted the nCSF model to the fMRI time series across the cortex, demonstrating its ability to capture a diverse range of responses. Figure 4 presents example fits from individual vertices in the V1 and TO regions (Fig. 4A), highlighting contrasting parameter values. The variance explained by the model for these locations was 74.72% for TO and 57.25% for V1, indicating that the nCSF model effectively captures the variance in the fMRI time series across both regions. Furthermore, the predicted fMRI time series closely align with the recorded data, as illustrated in Figure 4B. Notably, the nCSF curve in V1 exhibits a higher SF_p compared to the cortical location in TO (Fig. 4C).

Here, we report the parameters SF_p and normalized AUC as they are (1) they are the most robust to noise (see Section 3.1), (2) can be used as serve as a useful readout of the whole CSF (e.g., (Applegate et al., 1997)), and (3) can be consistently derived across a range of CSF parameterizations (Jigo et al., 2023).

The visual areas and eccentricity maps derived from pRF mapping for one participant are shown on the inflated cortical surface in the region near the occipital pole (Fig. 5, see Fig. S8 for the other participants). The eccentricity maps range from 0–5 degrees of visual angle, indicating the visual field representation. nCSF parameter and outcome estimates are projected on the inflated cortical surface; see Figure 5C-E.

The variance explained of the nCSF model is high across eccentricities and visual areas; see Figure 5C. The normalized AUC is high around the foveal region and generally decreases with eccentricity (Fig. 5D). SF_p is also higher near the foveal region compared to the parafoveal region and decreases as a function of eccentricity (see Figs. 5E and 6B), consistent with previous studies (Aghajari et al., 2020; Broderick et al., 2022; Henriksson et al., 2008). We observe that these systematic changes are similar between visual areas and participants (see Fig. S6).

To determine how nCSF parameters vary with eccentricity in V1, we used a linear regression analysis to fit the slope of eccentricity versus nCSF parameters. This analysis was performed separately per participant, only including vertices where both the nCSF and pRF model fits had a variance explained >30%. To correct for volume-to-surface upsampling, we determined the degrees of freedom used in calculating the t-statistics by taking the number of cortical locations (vertices) divided by the upsampling factor.

Overall, a decrease in SF_p and normalized AUC can be observed in the nCSF curves for different eccentricity bands (Fig. 6A for an example participant; Fig. S6 for all participants), displaying a shift leftward (toward lower spatial frequencies) as well as a general decrease in the area under the curve. SF_p decreases as a function of eccentricity in V1 for all participants, p < 0.001 (see Fig. 6B and Table S2). The normalized AUC also decreases with eccentricity in V1 for all participants, p < 0.001 (see Fig. 6C and Table S2).

We also compared how SF_P and normalized AUC differ from the horizontal and vertical meridian by performing an independent-sample t-test, separately per participant (again only including vertices where both the nCSF and pRF model fits had a variance explained >30% and correcting for volume-to-surface upsampling). Vertices were assigned to the horizontal or vertical meridian if they were within 15 (radial) degrees of the meridian, and had eccentricity values between 1 and 5 degrees (to avoid ambiguous polar angles around the fovea, following Himmelberg et al., 2022). Additionally, we found different nCSF curves for the horizontal and the vertical meridian (see Fig. 6D, example participant), and lower values for SF_P and normalized AUC in the vertical meridian compared to the horizontal meridian (p < 0.001 for all participants except P1, see Fig. 6E, F).

First, we compared nCSF curves across the cortical hierarchy. The nCSF curves are different across the cortical hierarchy, for example between V1 and TO. This is exemplified in Figure 7A, which shows the nCSF curves for one participant (for a comprehensive view of all participants and visual areas, see Fig. S6). Notably, the center of the nCSF curve for TO is further left (lower SF_p) than V1. Next, we compared SF_p and normalized AUC across eccentricity and the cortical hierarchy. Extrastriate regions tend to show a decrease in SF_p and normalized AUC with eccentricity, similar to V1 (see Fig. 7B, 7C for an example participant, and Fig. S7 for all other participants). However, this effect is less pronounced in extrastriate regions compared to V1. Additionally, there is an overall decrease across participants in SF_p and AUC moving up the cortical hierarchy (see Fig. 7E, 7F). Lastly, we investigated the goodness-of-fit of the nCSF model. Across all visual areas and participants, the nCSF model effectively captures the variance in the fMRI data (variance explained > 30%); see Figure 7D.

We introduce the neural contrast sensitivity function (nCSF) which is both a concept and a method. As a concept, the nCSF describes the sensitivity of neural populations as a function of spatial frequency and contrast. As a method, we model the nCSF in the human visual cortex using fMRI and using an approach similar to pRF modeling (Dumoulin & Wandell, 2008). The nCSF parameters were estimated from the fMRI data using a biologically-inspired, model-based approach similar to the pRF method. We model nCSF properties using a combination of the CSF (Chung & Legge, 2016) and CRF (Boynton et al., 1999).

We validated the nCSF method using simulations, showing that the key parameters of interest—spatial frequency peak and area under the curve—can be estimated accurately and are robust to noise. Importantly, the nCSF model effectively captured the variance in the fMRI data across the visual cortex. Consistent with prior fMRI studies (Aghajari et al., 2020; Goulet & Farivar, 2024; Henriksson et al., 2008), we found that the peak spatial frequency decreases as a function of eccentricity. In addition, total sensitivity (i.e., area under the curve) decreases with eccentricity. This effect was similar across the visual cortex, but more pronounced in the early visual cortex (V1, V2, V3) versus the late visual cortex (e.g. LO and TO). We also find that peak spatial frequency and area under the curve are lower in the vertical than horizontal meridian which is in line with findings in psychophysics (Himmelberg et al., 2020, 2022).

The nCSF concept combines notions of ophthalmology (Campbell & Gubisch, 1966; Le Grand, 1935; Michael et al., 2011), single neuron recordings (Albrecht & Hamilton, 1982; Levitt et al., 1994; Sclar et al., 1990), and neuroimaging (Aghajari et al., 2020; Boynton et al., 1999; Henriksson et al., 2008). In ophthalmology, the nCSF refers to processing in the retina and visual hierarchy but independent of optical aberrations. Though we measure brain responses directly, optical aberrations may still contribute. However, different nCSF properties in the same visual field location across the different visual field maps cannot be attributed to optical aberrations. Similarly, variations across eccentricity and polar angle are not likely due to optical aberrations, but due to known neural processing differences. Our nCSF concept also differs from single neuron recordings. For example, a population of neurons and the hemodynamic response properties contribute to the signal. These considerations are similar to the considerations underlying the concept of the pRF (Dumoulin & Knapen, 2018; Dumoulin & Wandell, 2008). We measure the nCSF using fMRI, but, in principle, the nCSF is independent of the measurement modality and can also be measured with other imaging techniques, such as invasive human electrodes (Harvey et al., 2013; Yoshor et al., 2007) or MEG (Eickhoff et al., 2024; Kupers et al., 2021).

Many factors may lead to changes in nCSF properties, for example using different stimuli or tasks. Different stimuli, for example different temporal characteristics, may elicit responses from different neural populations and therefore result in different nCSF. Here, we used gratings, but the same nCSF model could be fit and applied with any stimuli which covers the range of spatial frequencies and contrasts of interest; allowing comparisons across a large range of potential experiments. Likewise, attention influences the behavioral CSF, and therefore different tasks may influence the nCSF (Cameron et al., 2002; Jigo & Carrasco, 2020). Thus, the stimuli and tasks are likely to influence the neural population driving the responses underlying the nCSF model.

Furthermore, the nCSF model is based on a psychophysical model of the CSF, which can be parameterized using many different functions. We selected the asymmetric logarithmic parabolic function (aLP, Chung & Legge, 2016); fixing the width of the left side of the function, giving a 3 parameter model (CS_p, SF_p, width_r). We selected this model as a balance between flexibility and simplicity. The advantage of the aLP over symmetrical models is that it can capture the slower decrease in sensitivity to lower spatial frequencies compared with the more rapid decrease toward higher spatial frequencies, with the same number of parameters. Furthermore, disorders affecting the CSF typically (but not always, see below) affect it in characteristic ways, with either (1) a uniform drop in contrast sensitivity (translation down the y axis), (2) uniform drop in spatial frequency (translation leftward), or (3) selective sensitivity loss at high spatial frequencies (Chung & Legge, 2016). The aLP can appropriately capture all these changes, with a small number of parameters.

However, alternative CSF equations would likely produce similar results (Watson & Ahumada, 2005). Jigo et al. (2023) found that of the nine different functions, the “YQM” model (from Yang et al., 1995) provided the best fit to the data, but the aLP and several other models performed similarly. Previous fMRI studies of spatial frequency sensitivity have included a log-gaussian (Aghajari et al., 2020) and a symmetrical log parabolic function (Goulet & Farivar, 2024). Given the similarities between the functions (Watson & Ahumada, 2005), it is likely that these different models will perform similarly, especially on core parameters such as peak spatial frequency and total sensitivity (area under the curve).

It is likely the optimal form of the model will vary depending on use case. For some conditions, a two-parameter model may be appropriate (Chung & Legge, 2016). In some specific patient groups, a more sophisticated model will be necessary, for example selective loss of sensitivity at low spatial frequencies (cataracts) or selective loss at middle spatial frequencies (multiple sclerosis) would not be captured well by the three-parameter form of the aLP.

Just as there are alternative functions to describe the CSF component, there are alternative functions to describe the CRF. Here, we used a modified form of the Naka-Rushton function, with a single slope parameter allowing us to capture the nonlinear responses to increases in contrast. Many alternative functions have been used to fit the CRF, including the two-parameter form the Naka-Rushton function (Boynton et al., 1999) as well as linear, logarithmic, and exponential functions (e.g., Albrecht & Hamilton, 1982) and a binary step function (Goulet & Farivar, 2024). We selected the Naka-Rushton function as it is often used to model the CRF and is biologically plausible, approximating the gain control known to occur in cortical neurons (Heeger, 1992). For fMRI data, there is evidence that sustained adaptation promotes nonlinearities in the CRF, while unadapted responses will be more linear (Vinke et al., 2022). Hence, in some conditions, a linear model might be more appropriate, depending on the duration of the stimulus. Similarly, using stimuli which move more slowly through contrast space will likely result in more accurate estimates for the CRF (see CRF Simulation and Estimating the slope_CRF sections in Supplementary Materials).

In summary, our implementation of the nCSF offers a robust and accurate fit to the data while maintaining simplicity and interpretability. However, it is not the only possible implementation. In the same way that the pRF model (Dumoulin & Wandell, 2008) has been extended, we expect that the nCSF model will also be adapted and refined for different applications. Future models could benefit from integrating established computational motifs such as suppression (Zuiderbaan et al., 2012), response compression (Kay et al., 2013), and divisive normalization (Aqil et al., 2021) to enhance their predictive power. In studies involving patient populations where the shape of the CSF may vary or be less predictable, model-free approaches could be especially valuable for capturing unique characteristics and providing a more flexible framework (Carvalho et al., 2020; Greene et al., 2014; Lee et al., 2013).

The psychophysical CSF is widely used to assess perception in health and disease (Chung & Legge, 2016; Hou et al., 2010; Pelli & Bex, 2013). However, neural CSFs can vary across different cortical regions and levels of the processing hierarchy, and these variations may not always align with perceptual CSFs. This discrepancy highlights the need to investigate both perceptual and neural measures, as deficits might be present in one domain but not the other. The nCSF model provides a useful starting point for future research in this area, helping to better understand the diverse characteristics of neural computations across different contexts, populations, and conditions.

We introduce a novel method to measure the nCSF in the human visual cortex. We provide a quantitative, validated framework and show how nCSF properties vary across visual areas and with eccentricity. This method can be applied to both healthy and clinical conditions, and provide novel insights into the cortical organization underlying perception.

The code generated for this paper is available here: https://github.com/spinoza-centre/prfpy_csenf

The data are considered personal data pursuant to the General Data Protection Regulation (GDPR) and can only be shared based on and subject to the policies of the Royal Netherlands Academy of Arts and Sciences (KNAW). The datasets are available from the corresponding author on request. Considering the requirements imposed by law and the sensitive nature of personal data, any requests will be addressed on a case-by-case basis, subject to a data usage agreement.

S.O.D. and C.R. conceived the study. C.R. collected the data. C.R., M.D., and J.A.D. analyzed the data. C.R. wrote the paper. M.D., M.C.J., and S.O.D. edited the paper. S.O.D. and M.C.J. provided supervision.

The authors declare no competing interest.

This research was supported by Dutch Research Council Vici Grant 016.Vici.185.050 (to S.O.D.).

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00469.