Pupil-Linked Arousal Modulates Precision of Stimulus Representation in Cortex

Participants

Thirty-two healthy adult volunteers (20 female, 12 male, age: 19–31 years) with normal or corrected-to-normal vision participated in this study, which was approved by the local medical ethics review committee (CMO Arnhem-Nijmegen, the Netherlands). All participants provided informed written consent and received monetary compensation for their participation. The sample size (N = 32) was based on a power calculation for detecting a reliable correlation between decoded uncertainty and behavioral variability using data from a previous study with a similar design (van Bergen et al., 2015; power = 0.8; α = 0.05). Participants were included based on their ability to perform the task, which was assessed in a separate behavioral training session prior to the experimental sessions.

Imaging data acquisition

The MRI data were collected using a Siemens 3 T MAGNETOM Prisma Fit scanner and a 32-channel head coil at the Donders Centre for Cognitive Neuroimaging in Nijmegen, the Netherlands. The data were analyzed previously for a different purpose (Geurts et al., 2022). Each scan session started with the collection of a T1-weighted image [3D MPRAGE; repetition time (TR), 2,300 ms; inversion time (TI), 1,100 ms; echo time (TE), 3 ms; flip angle, 8°; field of view (FOV), 256 mm × 256 mm; 192 sagittal slices; 1 mm isotropic voxels] and B0 field inhomogeneity maps (TR, 653 ms; TE, 4.92 ms; flip angle, 60°; FOV, 256 mm × 256 mm; 68 transversal slices; 2 mm isotropic voxels; interleaved slice acquisition). Functional MRI data were acquired using a multiband accelerated gradient-echo EPI sequence, with 68 transversal slices covering the whole brain (TR, 1,500 ms; TE, 38.60 ms; flip angle, 75°; FOV, 210 mm × 210 mm; 2 mm isotropic voxels; multiband acceleration factor, 4; interleaved slice acquisition).

Pupil data acquisition

Pupillometry data were acquired using an SR Research EyeLink 1000 system. Pupil size was sampled at 1 kHz. Pupil recordings were collected for 62 out of 64 sessions and only partially (4–12 runs out of a total of 10–13) for 11 of these sessions, due to technical difficulties.

Experimental design

Participants performed an orientation estimation task (Fig. 1) inside the MRI scanner. They were instructed to maintain fixation on a black-and-white bullseye target (radius, 0.375°) presented at the center of the screen throughout each run. Runs consisted of 20 trials each (trial duration, 16.5 s; intertrial interval, 1.5 s) and started and ended with a fixation period (duration, 4.5 and 16 s, respectively). Each trial began with the presentation of an orientation stimulus (duration, 1.5 s), followed by a 6 s retention interval, and two 4.5 s response windows in which observers were prompted to report the orientation of the viewed grating and indicate their level of confidence in this orientation response. The stimuli were counterphasing sinusoidal gratings (contrast, 10%; spatial frequency, one cycle per degree; randomized spatial phase; 2 Hz sinusoidal contrast modulation), presented inside an annulus around fixation (inner radius, 1.5°; outer radius, 7.5°; contrast linearly decreasing over the inner and outer 0.5° of the annulus). Stimulus orientations were drawn pseudorandomly from a uniform distribution (0–179°) to ensure an approximately even sampling of orientation within any given run. During the first response window, participants reported the orientation of the viewed grating by rotating a black bar (length, 2.8°; width, 0.1°; contrast, 40%; initial orientation randomized across trials) presented at the center of the screen. During the second response window, participants indicated their confidence in this orientation judgment by sliding a white dot over a circular scale. The scale was a black bar of increasing width (contrast, 40%; bar width, 0.1–0.5°, linearly increasing) that was wrapped around fixation (radius, 1.4°). The mapping between confidence level and scale width (i.e., whether the narrow end of the scale indicated high or low confidence) was counterbalanced across participants, and the orientation and direction of the scale (i.e., whether the width increased in clockwise or counterclockwise direction), as well as the dot's starting position, were randomized across trials. For both response windows, the response bar (or scale) faded linearly over the last 1 s of the response window to indicate the approaching end of this window, and participants responded using two buttons (one for clockwise and one for counterclockwise rotation) on an MRI-compatible button box. Each trial was preceded by the fixation bullseye briefly turning black (duration, 0.1 s; timing, −0.5 s relative to stimulus onset) as a cue to stimulus onset. Participants performed a total of 22–26 task runs inside the scanner, divided over two sessions on separate days, and extensively practiced the task in separate behavioral sessions prior to the experiment (2–4 h in total).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Overview of the orientation estimation task. Participants were required to fixate the bullseye target in the center of the screen throughout each run. Each trial started with this bullseye briefly turning black, as a cue to stimulus onset. The visual stimulus was a counterphasing sinusoidal grating, presented in an annulus around fixation, and was presented for 1.5 s. After a 6 s delay (fixation interval), a black bar appeared in the center of the screen, and participants were required to report their orientation estimate by rotating this bar. Next, they reported their level of confidence about this estimate by sliding a dot over the (circular) confidence scale. Both response intervals lasted for 4.5 s, and the bar or scale started fading after 3.5 s to indicate the approaching end of the response window. A 1.5 s intertrial interval separated the response window from the next stimulus. Participants completed 20 trials per run.

Each scan session also included one or two functional localizer runs, in which flickering checkerboard stimuli (contrast, 100%; flicker frequency, 10 Hz; check size, 0.5°) were presented within the same annulus as the orientation stimuli. Checkerboard stimuli were presented in seven 12 s blocks interleaved with fixation blocks of equal duration. In a separate scan session, retinotopic maps of the visual cortex were acquired using standard retinotopic mapping procedures (Sereno et al., 1995; Deyoe et al., 1996; Engel et al., 1997).

Visual stimuli were generated by a MacBook Pro computer using MATLAB and the Psychophysics Toolbox (Brainard, 1997; Kleiner et al., 2007) and displayed via a luminance-calibrated projector (EIKI LC-XL100; screen resolution, 1,024 × 768 pixels; refresh rate, 60 Hz) on a rear-projection screen, which the participants viewed via a mirror mounted on the head coil.

Preprocessing of MRI data

Preprocessing procedures for functional imaging data are also described in Geurts et al. (2022) and reproduced here for convenience. Motion correction was performed with respect to the middle volume of the middle run of each session (the motion correction template) with FSL's MCFLIRT (Jenkinson et al., 2002). The motion-correction template was corrected for distortions due to B0 field inhomogeneities using the acquired field maps and registered to a high-resolution anatomical (T1-weighted) image acquired in the same session using epi_reg within FSL's FLIRT (Jenkinson and Smith, 2001). For coregistration of data across sessions, we created participant-specific anatomical templates by combining the anatomical reference images from the two separate sessions using FreeSurfer's mri_robust_template (Reuter et al., 2012), to which the single-session anatomical images were registered. All transformations were then combined and applied to the raw data. To remove slow drifts in the MRI signal, the transformed data were temporally filtered using FSL's nonlinear high-pass filter (Jenkinson et al., 2012) with a sigma of 24 TRs (two trials), which corresponds to a cutoff of ∼83 s. Residual motion effects were removed from the data through linear regression, using a set of 24 motion regressors derived from the motion parameters estimated by MCFLIRT.

The region of interest (ROI) for decoding (bilateral V1, V2, and V3, combined) was identified on the reconstructed cortical surface, obtained with FreeSurfer's cortical reconstruction algorithm (Dale et al., 1999), using single-participant retinotopic maps (see above, Experimental design). For further analysis, we selected the 2,000 voxels within this ROI that were most strongly activated by the functional localizer stimulus while surviving a lenient statistical threshold of p < 0.01, uncorrected. Voxel selection was performed for each participant individually in native space. Each voxel's time series were z-normalized with respect to corresponding trial time points in the same run. Finally, we obtained single-trial activation patterns by adding a 4.5 s temporal shift (to account for the hemodynamic delay) and then averaging over the first 3 s of each trial. Importantly, this time window excludes activity from the behavioral response window.

Decoding algorithm

To quantify the trial-by-trial precision of stimulus representations in the visual cortex, we applied a generative model-based probabilistic decoding algorithm (van Bergen et al., 2015; van Bergen and Jehee, 2018) to our data (see above, Preprocessing of MRI data, for voxel selection criteria). We provide a concise description of the decoding algorithm here and refer the interested reader to previous publications for further detail (van Bergen et al., 2015; van Bergen and Jehee, 2018). The decoding model describes the across-trial distribution of activation patterns as a multivariate normal, centered around a stimulus-specific mean that describes the tuning function of each voxel. Tuning functions were modeled as a linear combination Wf(s) of eight bell-shaped basis functions f(s)=[f1(s),…,f8(s)]T , each centered on a different orientation (Brouwer and Heeger, 2009). The basis functions are described as follows:fk(s)=max(0,cos(2πs−ϕk180))5, in which s is the orientation of the presented stimulus and ϕk is the center of the kth basis function. Basis functions were spaced equally across the full orientation space (0–179°) with the first one centered at zero. The basis functions are weighted by coefficients W, with Wik the contribution of the kth basis function to the tuning function of the ith voxel.

The variance around the stimulus-dependent mean activation pattern (i.e., the multivariate tuning function) is modeled by the covariance matrix:Ω=ρττT+(1−ρ)diag(τ2)+σ2WWT. The first component of the covariance matrix ρττT models global fluctuations shared between all voxels in the ROI. The second component(1−ρ)diag(τ2) , with τ2=[τi2]T , describes independent, voxel-specific variability (with variance τi2 for voxel i). The relative contributions of these two components are given by ρ. The third component σ2WWT reflects variance (with magnitude σ2 ) shared between voxels with similar orientation preference (given by WWT ).

The voxel tuning functions and covariance matrix together model the generative distribution of activation patterns:p(b|s,θ)=N(Wf(s),Ω), which is a multivariate normal distribution with mean Wf(s) and covariance Ω. θ={W,ρ,τ,σ} are the generative model's parameters.

For model training and testing, a leave-one-run-out cross-validation procedure was used to prevent double-dipping. The model's parameters were estimated on a dataset consisting of data from all but one task run. The trained model was then tested on the held-out run, and this procedure was repeated until all runs had served as a test run exactly once. The parameters were estimated in two steps: the coefficients W were first estimated by ordinary least squares regression, and then the covariance parameters (ρ,τ,σ) were estimated through numerical likelihood maximization (see van Bergen et al., 2015 for further details regarding model estimation procedures).

Model testing (“decoding”) consisted of calculating a posterior distribution p(s|b;θ^) over stimulus orientation s, conditioned on the response pattern b and estimated parameters θ^ . The posterior distribution is given by Bayes’ rule:p(s|b;θ^)=p(b|s;θ^)p(s)∫p(b|s;θ^)p(s)ds.

The stimulus prior p(s) was flat, reflecting the uniform stimulus distribution used in the experiment. The normalization constant ∫p(b|s;θ^)p(s)ds was computed numerically. The circular mean of the decoded distribution was taken as the decoder's estimate of the presented orientation (“decoded orientation”) and the squared circular standard deviation quantified the uncertainty in this estimate (“decoded uncertainty”).

Preprocessing of pupil data

Blinks and saccades were identified using the EyeLink software. Data recorded during saccades or <250 ms before (after) blink onset (offset) were removed. Missing or removed data were linearly interpolated. Data interpolated over >1,000 ms were removed at the end of the preprocessing procedure. If >50% of the data in a given trial were missing, the trial was excluded from pupil data analyses.

The pupil's responses to blinks and saccades were estimated and regressed out using a deconvolution approach developed by Knapen et al. (2016) and implemented by Urai et al. (2017). Specifically, the shape of blink- and saccade-triggered pupil responses was first estimated by fitting a finite impulse response (FIR) model to the data of each participant. The estimated response was then used to create a regressor in a linear regression analysis, and blink- and saccade-related effects were removed (separately for each run). Data were subsequently low-pass filtered using a third-order Butterworth filter with a cutoff of 4 Hz and downsampled to 100 Hz. Global effects in the data (cf. Knapen et al., 2016) were removed by fitting an exponential function to each run and using the residuals of this fit in subsequent analyses. Finally, pupil size was z-normalized per session to correct for differences in camera and lighting position between sessions.

We considered both the pupil size time series and its first derivative as indices of arousal. To estimate the derivative, we used a moving window with a width of 500 or 1,000 ms. Within this window, we fitted a linear function to the pupil time series. We took the slope of this fitted function as a measure of the rate of change in pupil size at the time point on which the moving window was centered.

Preprocessing of behavioral data

The error in the observer's orientation response on a given trial was calculated as the acute-angle difference between the presented and reported orientation. Participants generally performed well on the task, with a mean absolute orientation response error of 5.81 ± 1.29° (mean ± SD across observers). To correct for orientation-dependent biases in the response (shifts in mean response), two fourth-degree polynomials modeling response error as a function of stimulus orientation were fit to each participant's data (van Bergen et al., 2015; Geurts et al., 2022). The first polynomial was fit to trials with presented stimulus orientation between 0 and 89° and the second to trials with the presented stimulus between 90 and 179°. The residuals of these fits (“bias-corrected behavioral responses”) were used in subsequent analyses. Trials on which the bias-corrected response was more than three standard deviations away from zero were marked as guesses and excluded from all further analyses (0.91 ± 0.31 percent of all trials; mean ± SD across observers). Confidence ratings were z-scored per session to correct for potential between-participant or between-session differences in usage of the confidence scale. We excluded trials on which observers did not finish adjusting their orientation and/or confidence response before the end of the response window (2.75 ± 2.40 percent of all trials; mean ± SD across observers).

Statistical procedures

To benchmark our decoding approach, we quantified orientation decoding performance by calculating the circular equivalent of Pearson's correlation coefficient between the presented and decoded orientation across trials and for each individual observer. To quantify the effect at the group level, the single-observer correlation coefficients were Fisher-transformed, and a weighted average was computed. The weight for the correlation coefficient of observer i was calculated as wi=1vi , where vi corresponds to the variance of the Fisher-transformed correlation coefficient (Hedges and Olkin, 1985). This variance is given by vi=1ni−3 , in which ni represents the number of trials. Statistical significance of the group-averaged correlation coefficients was assessed using a z test. The results from this analysis, as well as the correlation between decoded uncertainty and distance to the nearest cardinal axis and behavioral variability (see below), were also reported as benchmarks in a previous study [Geurts et al. (2022), their Extended Data Fig. 2].

To assess whether decoded uncertainty predicts behavioral variability, trials were divided into 10 bins of increasing uncertainty for each individual observer. Mean decoded uncertainty was computed across all trials in each bin, and behavioral variability was quantified as the squared circular standard deviation of the (bias-corrected) behavioral errors in the bin. Multiple linear regression was performed to calculate the partial correlation coefficient for the relationship between decoded uncertainty and behavioral variability at the group level, with separate intercepts for each observer. Statistical significance was assessed by means of a t test. We performed a number of control analyses, in which we varied several analysis parameters. First, because there is no principled way for determining the number of bins to use in this analysis, we ran two additional analyses using 5 or 15 (instead of 10) uncertainty bins per participant. Second, because the strength of the link between decoded uncertainty and behavioral variability could vary across individuals, we performed a linear mixed-effects analysis (using MATLAB's lmefit function) in which both the intercept and the slope were modeled as random effects. This is in contrast to the multiple linear regression approach described above, which assumes that the intercepts are random variables while the slope is fixed. Third, to assess the influence of extreme values, we performed the analysis on different subsets of the data. The first subset excluded all trials for which the standard deviation of the decoded distribution was larger than 45°. In the second subset, we simply excluded from our analyses the observer who gave rise to the data point in the top-right corner of Figure 2B. Statistical significance was assessed by means of t tests on the partial correlation coefficients (r) or estimated slope (β) for the multiple regression and linear mixed effects analyses, respectively. In the linear mixed effects analysis, Satterthwaite approximation was used to estimate the effective degrees of freedom.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Overview of behavioral, physiological, and neural measures and behavioral benchmarking. A, The imprecision in the cortical stimulus representation was quantified using a probabilistic decoding technique. This decoding algorithm computes a probability distribution over stimulus orientation from single-trial activity patterns. The width of the distribution was taken as a measure of the degree of uncertainty in the cortical stimulus representation. Pupil size and the slope (the first temporal derivative) of the pupil signal were used to measure arousal over the course of each trial. Behavioral imprecision was measured as the variability in the observer's orientation judgments across trials. B, Decoded uncertainty predicts behavioral variability on a trial-by-trial basis. Each observer's trials were split into 10 bins of increasing decoded uncertainty. Behavioral variability was computed as the squared circular standard deviation of behavioral response errors in each bin. The partial correlation coefficient between decoded uncertainty and behavioral variability (controlling for differences in intercept between participants) was computed and found to be significantly positive (t₍₂₈₇₎ = 2.30; p = 0.011; r = 0.13; 95% CI = [0.019, 0.25]). Thus, when the cortical stimulus representation was more uncertain, behavioral imprecision was larger, as well. Note that the data is centered around zero because this is a partial correlation plot so interindividual differences in the mean have been removed. The fMRI results in B were also reported in Geurts et al. (2022), their Extended Data Figure 2.

The relationship between pupil size and decoded uncertainty or reported confidence was tested in two different ways. For the first set of analyses, we divided trials into three bins per observer, based on the level of decoded uncertainty (or reported confidence), and compared pupil size between the highest and lowest bin. We computed t values to quantify the difference between these bins. t values were computed for each observer individually and then averaged across observers. The analysis was performed for each time point within the window of interest; that is, from 1.5 s before stimulus onset to 1.5 s after stimulus offset. To assess statistical significance, threshold-free cluster enhancement (TFCE; Smith and Nichols, 2009) and permutation testing (1,000 permutations) were used. The familywise error rate (FWER) was controlled by comparing the true single-time point TFCE scores against the null distribution of the maximum TFCE score across time obtained through data permutation (Nichols and Hayasaka, 2003).

In the second set of analyses, we quantified the relationship between pupil size (or slope) and decoded uncertainty (or reported confidence) on a trial-by-trial basis. To do so, we computed Spearman's correlation coefficient for each data point in the pupil time series, from 1.5 s before stimulus onset to 1.5 s after stimulus offset. Spearman's correlation coefficient was used because there is no a priori reason to assume that the relationship should be linear. Correlation coefficients were computed for each observer individually and group-level correlation coefficients were calculated following similar procedures as for orientation decoding performance. Specifically, individual correlation coefficients were Fisher transformed, and a weighted average was computed with weights wi=1vi (Hedges and Olkin, 1985). For Spearman's correlation coefficient, vi is given by vi=1.06ni−3 , with ni representing the number of trials for observer i (Fieller and Pearson, 1961). As before, statistical significance was assessed using TFCE and permutation testing (1,000 permutations), and the FWER was controlled by comparing against the null distribution of the maximum TFCE value across time points.

For the visualizations in Figure 4D, pupil size was averaged over the stimulus presentation window on a trial-by-trial basis. Spearman correlation coefficients between decoded uncertainty and (mean) pupil size were computed per observer, Fisher-transformed and averaged as described above. z tests were used to assess significance both at the group level and for the example observer.

The relative effect size of arousal state on changes in uncertainty was determined as follows. Because the effect size cannot be determined directly from the observed relationship between arousal and uncertainty (decoded uncertainty reflects not only neural but also many non-neural sources of noise, including from the fMRI scanner), we instead relied on an indirect approach and compared the effect of pupil-linked arousal with that of a manipulation of stimulus orientation to acquire an understanding of its relative contribution. The impact of stimulus orientation on decoded uncertainty was quantified by computing the Spearman correlation coefficient between decoded uncertainty and the distance between the presented stimulus orientation and the nearest cardinal axis [cf. Geurts et al. (2022), their Extended Data Fig. 2B]. Correlation coefficients were computed per observer, Fisher transformed, and then averaged across observers as described above. The impact of pupil-linked arousal was summarized by first averaging pupil size over the stimulus presentation window on a per trial basis (cf. Fig. 4D) and then computing and averaging the Spearman correlation coefficients between these trial-by-trial values and decoded uncertainty as described in the previous paragraph. Effect sizes were defined as the absolute group-averaged Spearman correlation coefficient.

Code accessibility

All custom analysis code is available from the corresponding author upon request. Code for the probabilistic decoding algorithm can be found at https://github.com/jeheelab/.