Twelve healthy right-handed participants (five males and seven females, aged between 19 and 29 years) with normal or corrected-to-normal vision were recruited in the present study. The study was approved by the Institutional Review Board of the Department of Psychology and the Center for Biomedical Imaging Research at Tsinghua University. All participants provided written informed consent before the experiments.

Four circular flickering checkerboards with different contrast levels were presented to elicit BOLD responses. The visual stimuli were displayed via a mirror and screen system. The visual stimuli were presented on an LCD screen, which was viewed through an angled mirror attached above the head coil. As there was an upgrade of the visual stimulus system during the research, an Invivo system (Gainesville, FL, United States) was used for the first four participants (screen size = 64 cm × 40 cm, corresponding visual angle ≈ 19.7° × 12.4°, resolution = 1280 × 800, refresh rate = 60 Hz) and a Sinorad system (Shenzhen, Guangdong, China) was used for other participants (screen size = 89 cm × 50 cm, corresponding visual angle ≈ 24.6° × 14.0°, resolution = 1920 × 1080, refresh rate = 60 Hz). The visual stimuli consisted of four circular flickering checkerboards (radius ≈ 2.0°, checker size ≈ 0.3° × 0.3°, flickering rate = 6 Hz), which were presented in four respective quadrants of a uniformly gray background with their centers ∼4.4° away from the center of the screen. Since it is difficult to map out the visual area boundaries corresponding to the fovea (Zeki, 1969; Dougherty et al., 2003; Schira et al., 2009; Wandell and Winawer, 2011) and the meridians that are sensitive to eye movements (Engel et al., 1997; Noory et al., 2015), the visual field corresponding to fovea as well as the horizontal and vertical meridians were avoided to constrain the BOLD responses inside V1. The checkerboards with a contrast level of 25% (Michelson contrast, adaptation contrast) were presented during the adaptation period. After this period, the checkerboards with a contrast level of 100 or 6.25% (Michelson contrast, test contrast) were presented to elicit respective BOLD signal increments or decrements (Figure 1). Increments and decrements of the contrast level were introduced to a logarithmic scale (in octaves, one step equals 2-fold), and a 4-fold contrast level change was used to elicit the largest BOLD signal increments or decrements (Gardner et al., 2005b). In addition to the four circular checkerboards, a circular fixation dot (radius ≈ 0.1°) was placed at the center of the screen throughout the whole experiment.

Contrast ramps were introduced when the stimuli changed their contrast level (Tuan et al., 2008). To be precise, checkerboards with intermediate contrast levels were presented sequentially during the change, and the contrast level would increase or decrease linearly as time went by. There were seven intermediate states for each change, and each intermediate state took 1/60 s.

All participants were instructed to complete a localizer experiment and a contrast adaptation experiment in the present study. The localizer experiment was conducted to identify the activated voxels. The localizer experiment consisted of two types of blocks which were a 30-s rest block and a 30-s stimulus block. The two blocks were presented and alternated four times in each run. No visual stimulus was presented on the gray background during the rest block, and checkerboards with a contrast level of 100% were presented during the stimulus block.

The event-related contrast adaptation experiment was performed to examine the linearity of BOLD responses (Figure 1). The contrast adaptation experiment consisted of three runs, with each run lasting for 890 s. For each run, no stimulus was presented at the background during the initial 30 s, then the adaptation stimulus with a contrast level of 25% was presented for 60 s. After the adaptation phase, 24 trials were conducted sequentially with each trial including a test stimulus presentation followed by a 30-s adaptation stimulus presentation. The test stimulus was presented for 1, 3, or 6 s, with a contrast level of 100 or 6.25%. Four repetitions for each combination of stimulus duration and contrast level were included in each run, and the order of trials was randomized.

Throughout the whole experiment, participants were instructed to conduct a button-pressing task. Every 1–5 s, the white fixation dot that the participants fixated on turned red for 200 ms, and participants were required to report the change of color with a button press within 1 s.

Two different MRI systems were used during the experiment due to a mandatory upgrade of the system. A Philips Achieva 3 T MRI system was used for the first four participants and a Philips Ingenia CX 3 T MRI system was used for all other participants. Upon the MRI system upgrade, a serious artifact appeared in functional images using the new Philips Ingenia CX 3 T MRI system when the original Echo Planar Imaging (EPI) sequence that was used for the older Philips Achieva 3 T MRI system was applied to the new system. Therefore, a new EPI sequence with slightly different parameters was designed and applied under the guidance of professionals to adapt to the new Philips Ingenia CX 3 T MRI system. It should be noted that there were no systematic differences in the measurements between the new and old MRI systems.

The imaging data were recorded using a 32-channel radio-frequency coil. The whole-brain structural images were acquired using a T1-weighted Turbo Field Echo sequence (TR = 7.6 ms, TE = 3.7 ms, FOV = 23 × 23 cm², in-plane voxel size = 0.96 × 0.96 mm², in-plane matrix size = 240 × 240 pixels, 180 contiguous slices, slice thickness = 1 mm). Then the images of the slices that were perpendicular to the calcarine sulcus were collected, covering most areas of the occipital lobe extending from the occipital pole. The structural images for these slices were collected using a T1-weighted Turbo Field Echo sequence (TR = 2.2 s, TE = 13 ms, FOV = 23 × 23 cm², in-plane voxel size = 0.45 × 0.45 mm², in-plane matrix size = 512 × 512 pixels, 13–14 contiguous slices, slice thickness = 3 mm), and the functional images for these slices were acquired using an EPI sequence during the experiments (TR = 1 s, TE = 35 ms, flip angle = 90°, slice thickness = 3 mm). The EPI sequence that was used for the older Philips Achieva 3 T MRI system (FOV = 19.2 × 19.2 cm², in-plane voxel size = 3 × 3 mm², in-plane matrix size = 64 × 64 pixels, 13 contiguous slices, no Multi-band SENSE used) was slightly different from the EPI sequence that was used for the new Philips Ingenia CX 3 T MRI system (FOV = 22.0 × 22.0 cm², in-plane voxel size = 2.75 × 2.75 mm², in-plane matrix size = 80 × 80 pixels, 14 contiguous slices, Multi-band SENSE used, MB factor = 2). Lastly, additional empty scans were included at the beginning of each functional run, allowing the longitudinal magnetization to reach a steady state. The additional scans were excluded from the functional data that would be used for further analysis.

The functional images were head-motion corrected, followed by undergoing different pre-processing procedures for the localizer experiment and contrast adaptation experiment respectively. They were then registered to the structural images.

Firstly, head-motion correction was performed for the functional data via a two-step procedure: all EPI volumes were aligned to the first volume of the corresponding run and were then further aligned to the first volume of the time course of the localizer experiment. The 3dvolreg program in AFNI¹ (Cox, 1996) was used to conduct this procedure.

Secondly, the functional images for the localizer experiment were further pre-processed. Considering that the frequency range of BOLD signal fell between 0.0167 and 0.15 Hz (this frequency range was calculated based on the simulated BOLD time courses, see Supplementary Figure 1 for details), the linear trend for the BOLD time course was first removed, and a band-pass filter (high-pass frequency = 0.0125 Hz, low-pass frequency = 0.2 Hz) was then applied to remove the low and high-frequency noise. The 3dFourier program in AFNI was used for temporal filtering, removing the mean and the linear trend before each filtering.

Thirdly, the functional images for the contrast adaptation experiment were pre-processed by going through a series of procedures. As the initial period of the time course contained an ascending trend (the BOLD signal remained low in the first 30 s baseline, then increased to a higher level when the adaptation stimulus was presented, see Figure 2A), the linear trends of signal and noise were mixed together, and a complete removal of the linear trend would remove the signal that needed to be retained. Thus, the linear trends of the test period were removed after the first 90 volumes in each run were discarded. The simulation results suggested that the frequency range of the BOLD signal was approximately 0–0.15 Hz for the remaining time course in each run (Supplementary Figure 1), which was intermixed with the low-frequency noise. Instead of applying a band-pass filter on each run directly, a two-step filtering procedure was introduced. Specifically, the high-frequency noise was firstly removed with a low-pass filter (low-pass frequency = 0.2 Hz) in each run, followed by recombining the time courses of all trials of the same condition for three runs to form new time courses (Supplementary Figure 1). For instance, if a 1-s stimulus occurred 200 s after the start of a run, the trial for this stimulus that corresponded to volume 201–231 of this run would be selected to form the recombined time courses. The newly formed time courses with signals starting from the frequency that corresponded to the 12^th FFT component (12 cycles per time course) were then filtered with a high-pass filter, in which the high-pass frequency that equaled the frequency of the 10^th FFT component was applied to the newly formed time courses to remove the low-frequency noise. The frequency of the 10^th FFT component was calculated according to the time length of the newly formed time course. Each recombined time course contained the time courses of 12 trials, and a total of 6 recombined time courses were formed for the 6 conditions (2 BOLD response type × 3 stimulus duration), respectively. As each trial contained a test stimulus presentation followed by a 30-s adaptation stimulus presentation, the time length for the newly formed time courses under the 1-s stimulus condition was 31 × 12 = 372 s, and the frequency of the 10^th FFT component was 10 cycle / 372 s = 0.0269 Hz, which would be applied to remove the low-frequency noise. After the temporal filtering procedure, the measured BOLD responses of the six conditions were calculated for every voxel by averaging the time courses of the 12 trials. To unify the time length of the measured responses across all conditions, the first 25 volumes (corresponding to 25 s) of the measured responses were used in the linearity analysis, signal-to-noise ratio (SNR) analysis, and for calculating r² of the contrast adaptation experiment (see “Linearity Analysis,” “Signal-to-noise Ratio (SNR) Analysis,” and “Activated Voxels Selection” sections).

Finally, the functional images were registered to the high-resolution structural images using the mrAlign tool in mrTools (Gardner et al., 2018). The first volume for the functional images of the localizer experiment was registered to the structural images for the corresponding slices, where these structural images were further registered to the whole brain structural images.

In addition, all other pre-processing procedures were conducted by custom programs written in MATLAB R2016b (Mathworks, Natick, MA, United States). The programs (i.e., cbiReadNifti.m and cbiWriteNifti.m) in mrTools (Gardner et al., 2018) were used to read and write NIfTI files.

To choose voxels that were significantly activated by the visual stimuli within V1, the activated voxels were selected following three steps for each participant (Gardner et al., 2005b): (1) the boundaries of V1 were drawn based on a separate retinotopy experiment; (2) the voxels that were significantly activated by the visual stimuli were determined using the BOLD time courses of a block design localizer experiment (see “Experimental Design” section for the procedure of localizer experiment); (3) combining the results of the previous two steps, clusters of significantly activated voxels were chosen within V1. The result of the contrast adaptation experiment was not involved in the selection of activated voxels, avoiding the interaction between activated voxel selection and linearity analysis.

The boundaries of V1 in each participant were determined in a separate retinotopy experiment using the traveling-wave method (Engel et al., 1994, 1997; Wandell et al., 2007). Six runs of rotating wedge stimuli were presented to map retinotopic organization with respect to visual polar angle, and four runs of moving ring stimuli were presented to map retinotopic organization with respect to visual eccentricity. Each run took 252 s and contained 10.5 cycles. After removing the imaging data for the first 0.5 cycles of each run, the BOLD time courses for the runs with the wedge stimuli and ring stimuli presented were averaged separately. Subsequently, the averaged time courses in each voxel were fitted to the sinusoids at the stimulus frequency. A high correlation between the actual and fitted time courses suggested that the voxel was more likely to be within the visual cortex. In addition, the fitted sinusoids for some voxels with high correlation reached their peak value when the visual stimulus appeared at vertical meridians, and the boundaries of V1 were drawn on these voxels on the flattened cortical surfaces. The data analysis was conducted with mrTools, and the detailed data analysis procedures are available online². The flattened cortical surfaces were created by mrTools based on the segmentation done by FreeSurfer³ (Fischl, 2012).

To determine significantly activated voxels, an index (r²) was calculated for each voxel’s time course that was acquired from the localizer experiment, indicating what percentage of the variance in the time course was accounted for by the visual stimuli (Gardner et al., 2005b). To calculate r², an HRF was first estimated for each voxel’s time course using a deconvolution method (also called the finite impulse response model) (Dale, 1999; Hinrichs et al., 2000). Based on our previous study (Gardner et al., 2005b), the deconvolution method was chosen instead of the generally used GLM method as it does not presume any particular shape for the time course of the HRF. The goodness-of-fit (r² in this study) was then estimated, which is equal to the percentage of variance in the original time course that is accounted for by the estimated time course. The estimated time course was the convolution of the estimated HRF and the time course of the visual stimuli. The value of r² ranges from 0 to 1. The higher the value of r², the more variance was accounted for by the estimated time course. The value of r² is equal to 0 if the estimated time course does not account for any of the variance, and the value of r² is equal to 1 if the estimated time course accounts for all the variance. The calculation of r² was implemented by mrTools⁴. To determine the significance level of a certain r² value, the r² distributions of randomized time courses were constructed. For each participant, the BOLD time courses for all voxels in retinotopically defined V1 were extracted after the head-motion-correction procedure. This was followed by random shuffling of the time points in each BOLD time course, creating 100000 randomized time courses in total. After these randomized time courses went through the pre-processing procedures (i.e., removing the linear trend, filtering), the values of r² were computed for each time course, and these values were used to form the randomized distribution of r². The false positive rate for the value of r² for each voxel was determined based on the randomized distribution. For example, if the value of r² was ranked as the 99.9% highest in the randomized distribution, then its false positive rate would be equal to 0.001. The voxels that had a false positive rate below 0.001 were shown for one of the participants in Supplementary Figure 2. The false positive rates for the voxels in V1 were corrected using the false discovery rate (FDR) approach (Benjamini and Hochberg, 1995). A voxel was considered significantly activated when its corrected false positive rate was smaller than 0.001. Another randomized distribution using simulated time courses with temporal autocorrelation was generated to verify the reliability of this result (see Supplementary Material). When computing the randomized distributions, programs (i.e., cbiReadNifti.m and cbiWriteNifti.m) in mrTools were applied to read and write NIfTI files, while all other processes were conducted by custom programs written in MATLAB R2016b.

Clusters of significantly activated voxels were selected within V1 of each participant. The size of clusters was defined based on the size of the visual stimuli according to previous studies (Horton and Hoyt, 1991; Sereno et al., 1995). In addition, the visual stimuli in the localizer experiment elicited a BOLD signal increment in all selected voxels. Furthermore, the r² values for every selected voxel were calculated using the time courses of the contrast adaptation experiment to verify that the selected voxels were responding to the visual stimuli in the contrast adaptation experiment. The result showed that 98.7% of the selected voxels had their r² value above the chance level (false positive rate = 0.001), suggesting that most selected voxels responded to the visual stimuli in the contrast adaptation experiment.

To obtain a voxel-wise test of linearity, the similarity between measured BOLD responses and their linear predictions was assessed in every selected voxel by calculating the value of a similarity index, the Dice index (Dice, 1945; Cha, 2007). The linear predictions for the measured BOLD responses were calculated by adding the BOLD responses to the shorter stimuli with their temporally shifted copies (Boynton et al., 1996). For instance, to predict the BOLD signal increment to the 6-s stimulus, the BOLD signal increment to the 3-s stimulus was summed up with the copy of its shifted self by 3 s. In the present study, the linear predictions of responses to the 3-s stimulus were made by responses to the 1-s stimulus, and the linear predictions of responses to the 6-s stimulus were made by responses to the 1-s stimulus and 3-s stimulus, respectively. It is important to underline that the linear prediction should approximate measured BOLD response when the linearity holds, whereas deviations should occur when the BOLD responses behave in a non-linear manner. Therefore, the more the BOLD response behaved in a linear way, the higher the degree of similarity between the measured BOLD response and its linear prediction. Previous studies suggested that the difference in the amplitude of BOLD time courses was the major indicator of non-linearity (Birn et al., 2001; Soltysik et al., 2004). Hence, the Dice index was chosen as the suitable indicator as it is sensitive to the difference in the amplitude between time courses. Other similarity indices (e.g., the correlation coefficient) tend to describe the similarity in the direction of vectors, thus they can only describe the similarity in the shape but not the amplitude between time courses, making them unsuitable for the present study. Assuming a measured BOLD response (M) and its linear prediction (P) with n time points, their Dice index was calculated using the formula below (the subscript t refers to the t^th time point):

The value of the Dice index ranges from −1 to 1: the value is equal to 1 if and only if P = M, and the value is equal to −1 if and only if P = −M. The value increases if the similarity between P and M increases. Wilcoxon signed ranks tests were performed to assess whether there was a difference in the value of the Dice index across the BOLD response types and stimulus durations. In these tests, one data point referred to the Dice index of one selected voxel. These tests were implemented by the “signrank” function in MATLAB R2016b. To assess how many voxels behaved in a linear manner, the percentages of voxels that had their Dice index above chance level (false positive rate = 0.05) were calculated under each condition. The Dice thresholds for chance level were calculated based on the randomized distributions that were established by computing the Dice index for randomly shuffled time courses (see “Activated Voxels Selection” section for randomly shuffled time courses).

The voxel-wise overestimation and underestimation patterns were examined by estimating the contrast index for the amplitude of the HRFs in each selected voxel (Birn et al., 2001; Soltysik et al., 2004). HRFs were fitted to the measured BOLD responses under each condition, then the contrast index for the amplitude of the HRFs was calculated across stimulus durations of the BOLD signal increments and decrements, respectively. Firstly, a two-gamma HRF was fitted to the measured BOLD responses under each of the six conditions (see “Data Pre-processing” section for measured BOLD responses). The measured BOLD responses were either increment (for the 100% contrast condition) or decrement (for the 6.25% contrast condition), thus the HRF was either a positive HRF or a negative HRF. The positive HRF was modeled by a canonical HRF, and the negative HRF was modeled by zero minus canonical HRF. The formula for canonical HRF is given below (Lindquist et al., 2009):

The two α parameters were replaced by μβ + 1 in practice, where μ is the time each Gamma function reaches its peak. The formula has seven free parameters with ranges restricted to obtain a reasonable fit: μ₁ (from 0.5 to 8), μ₂ (from 4 to 16), β₁ (from 0 to 5), β₂ (from 0 to 5), c (from 0 to 1), t_onset (from 0 to 5), and A (from 0 to 20). After convolving the HRF with the time course of the visual stimuli (which was modeled as a boxcar function), the time course was fitted to the measured response using the trust-region-reflective algorithm (via “lsqcurvefit” function in MATLAB R2016b). An HRF with more flexible parameters was also applied to prove that the deviation patterns did not result from the restricted parameters. The ranges of parameters for the flexible HRF are μ₁ (from 0 to 24), μ₂ (from 0 to 24), β₁ (from 0 to positive infinity), β₂ (from 0 to positive infinity), c (from 0 to 1), t_onset (from 0 to 24), and A (from 0 to positive infinity). Secondly, the amplitude of HRF was estimated as the maximum of positive HRF or the minimum of negative HRF (note that both the maximum and minimum represent the peak of HRFs, thus the undershoot or overshoot after the peak was not considered). Thirdly, the contrast index for the amplitude of HRFs across stimulus durations was calculated for the BOLD signal increments and decrements separately. The formula for the contrast index is given below:

The amplitude of the HRFs for the short stimulus (Amp_short) was compared with the amplitude of HRFs for the longer stimulus (Amp_long). To be precise, the amplitude of the HRFs for the 3-s stimulus condition was compared with the amplitude of the HRFs for the 6-s stimulus condition, and the amplitude of the HRFs for the 1-s stimulus condition was compared with the amplitude of the HRFs for the 3-s and 6-s condition. The HRF for the short stimulus corresponded to the linear prediction, whereas HRF for the longer stimulus corresponded to the measured BOLD responses. Note that in the event when the BOLD response behaved in a linear manner, the linear prediction should approximate the measured BOLD responses, which yielded similar HRFs, and the contrast index should be close to 0. On the contrary, when an overestimation (or underestimation) pattern occurred, the contrast index should be above (or below) 0. The percentage of voxels that had a contrast index above 0 was calculated under each condition. For each distribution of the contrast index, a Kolmogorov–Smirnov one-sample test was performed to assess if the distribution significantly deviated from a normal distribution with its mean equaled to 0 and with the same standard deviation as the distribution for the contrast index. This test was implemented by the “kstest” function in MATLAB R2016b, and one data point in this test referred to the contrast index of one selected voxel. Like previous studies that used the ratio for the amplitudes of the BOLD responses to assess the linearity (Birn et al., 2001; Birn and Bandettini, 2005), the HRF amplitude ratios across stimulus durations were also calculated for the BOLD signal increments and decrements separately. To calculate the HRF amplitude ratio, the amplitude of the HRFs for the short stimulus was divided by the amplitude of HRFs for the longer stimulus. It should be noted that the distribution of the HRF amplitude ratio is asymmetric and unbounded for positive values, while the distribution of the contrast index is symmetric and ranges from −1 to 1, making the contrast index more appropriate for the statistical analysis. Therefore, the contrast index was used as the main indicator of the deviation pattern in the present study.

To support the validity of linearity analysis, the deviation pattern across the high SNR group and low SNR group were compared for the BOLD signal increments and decrements separately. For each selected voxel, the pre-processed time courses of all trials for the BOLD signal increments and decrements of three runs were recombined separately, and the SNR was calculated for each of the recombined time courses. The time course of each trial was unified to the first 25 volumes (seconds) starting from stimulus presentation, which yielded a task frequency of 25 s per cycle. The SNR was calculated by dividing the magnitude of the BOLD signal at task frequency (0.04 Hz) to the mean magnitude of the BOLD signal at 0.33–0.50 Hz, which represented the frequency range for high-frequency noise (Sun et al., 2013; Kuriki et al., 2015). A Wilcoxon signed ranks test was performed to compare SNR across BOLD signal increments and decrements. To explore the effect of SNR, the selected voxels were divided into a low and a high SNR group by the median SNR (two different divisions were created for the BOLD signal increments and decrements separately). Kolmogorov–Smirnov tests were performed to assess if the distributions of the contrast index differed across high SNR and low SNR groups using the “kstest2” function in MATLAB R2016b. In these tests, one data point referred to the contrast index of one selected voxel. Other noise bands (e.g., all frequencies except the task frequency) were also applied to prove that the choice of noise bands did not affect the deviation patterns in two SNR groups.

Limited by the length of the experiment, only three contrast levels (100, 25, and 6.25%) were chosen for the present study. The three contrast levels were chosen based on a logarithmic scale (Ohzawa et al., 1985), maximizing the difference between the contrast levels. Future studies are advised to apply more contrast levels in the experiment which can help to better understand the relationship between the stimulus intensity and the linearity of BOLD responses.

While previous studies used the averaged time courses for linearity analysis (Birn and Bandettini, 2005; Tang et al., 2009), the present study used the time course of each selected voxel. Although the BOLD signal increment showed a consistent overestimation pattern, the BOLD signal decrement showed an underestimation-to-overestimation pattern across the time courses of the selected voxels, which was different from the underestimation that was revealed in the averaged time courses. As the SNR of BOLD signal decrement was lower than the SNR of BOLD signal increment, the underestimation-to-overestimation pattern in the BOLD signal decrement could potentially be a result of the low SNR. To control the effect of the SNR, the deviation patterns of the time courses with similar SNR distribution were computed, and the observed deviation patterns for both BOLD signal increments and decrements were similar to the pattern for all selected voxels. This result suggested that the low SNR in the BOLD signal decrement was less likely to be the cause of the variations in deviation pattern. Nevertheless, the effect of SNR could not be completely excluded, and obtaining data with a higher SNR could increase the reliability of the results.

Apart from the BOLD signal decrement in the present study, previous studies also demonstrated another type of stimulus-induced BOLD signal decrease that was relative to the baseline condition, which was named as the negative BOLD responses (Fransson et al., 1999; Shmuel et al., 2002, 2006; Boorman et al., 2010; Gonzalez-Castillo et al., 2012; Klingner et al., 2015). For example, visual stimuli could elicit negative BOLD responses around the classical positive BOLD responses in the human visual cortex (Shmuel et al., 2002). Testing the linearity presumptions under the experimental settings for negative BOLD responses will help to reach a more general conclusion about the linearity of signal decrements in the future. In addition, dividing the BOLD responses into BOLD signal increments and decrements only focuses on one aspect of the temporal dynamics of BOLD responses. The BOLD responses demonstrate widely diverging temporal signatures across brain regions (Gonzalez-Castillo et al., 2012), thus further explorations are invited to examine the linearity presumption beyond the scope of BOLD signal increments and decrements.