Sixteen healthy subjects (8 males, 8 females; ages: 21-63, mean=41.2) completed the fMRI study after providing signed informed consent, approved by the Institutional Review Board at the National Institutes of Health (Combined Neurosciences White Panel). Participants were recruited and screened to exclude ferromagnetic implants, psychoactive medications and major medical problems, substance use disorders (other than nicotine), and neurological or psychiatric disorders (including eating disorders), as assessed by the Structured Clinical Interview for the Diagnostic and Statistical Manual of Mental Disorders (DSM-IV) 35. Data from two subjects were excluded from final analysis due to their failure to achieve criteria described below.
MRI data acquisition sequences
MRI data was collected on a 3T Magnetom Prisma scanner (Siemens Medical Solutions USA, Inc., Malvern, PA) using a 32-channel head coil. A standard echo planar imaging (EPI) sequence (sequential interleaved acquisition, repetition time 2.13s, echo time 30ms, flip angle α=79˚, 64×64 pixels in-plane resolution, 40 slices, slice thickness 3.5 mm, voxel dimensions 3×3×3.5 mm3, field of view 192×192 mm2) was used for functional scans. The three-dimensional, magnetization-prepared rapid gradient-echo 36 (3D MP-RAGE; TR/TE = 2400/2.24ms, FA = 8 deg, TI=1060ms), and Sampling Perfection with Application optimized Contrasts by using different flip angle Evolutions (SPACE, Siemens; TR/TE = 3200/564ms) pulse sequences were used to acquire high-resolution anatomical brain images with 0.8 mm isotropic voxels field-of-view (FOV) = 240 × 256 mm, matrix = 300 × 320, and 208 sagittal slices.
All participants underwent two fMRI scan sessions, each consisting of three runs (~13 minutes each) of a choice-reward task and two resting state runs (~10 minutes each, prior to and after the choice-reward fMRI task). Both choice-reward task and resting state experiments were developed using E-Prime software (Psychology Software Tools, Pittsburgh, PA, version 2.0.10 Professional). Stimuli were presented on a black background under dimmed room lighting using a liquid-crystal display screen (BOLDscreen 32, Cambridge Research Systems; UK). During the experiment the E-Prime computer communicated with the eye-tracker and the fMRI computers, either to send stimulus onset triggers, or to receive stimulus to start and stop the experiment. During the resting scan subjects were asked to keep their eyes open and to stare at a white cross on a black background and minimize eye movements. The choice-reward task is described in detail in Supplementary Methods. Briefly, participants were given 10$ worth of points to start the game. The reward value of the trial (low, 5 or 25 cents; high, 1 or 2 dollars, or Null: 0 money) was shown, and after three seconds, they had to make a choice between two shapes (a triangle or a square) presented side by side separated by a color (red or a green) circle at the center of the field of view. The green circle corresponded to a “GO” trial for which the participant had to make a choice between shapes by pressing the right or left button within a 2-seconds response window and their reaction times were measured to assess their association with blink latency. The red circle corresponded to a “NOGO” trial for which participants had been instructed not to respond. Four to six seconds after their response (or no response), a feedback message informed them of their choice outcome (win or lose money, or no gain/loss). For both, GO and NOGO trials winning was preset at 55% and loosing at 45%, independent of the participant’s choices, such that all participants would earn the same amount of points at the end of the task unless they missed a response or made a response to NOGO, situations where they were punished by subtracting 25 cents from their points.
Eye tracking data acquisition and analysis
During the fMRI scans, we measured pupil diameter and eye-blinks with an ASL (Applied Science Laboratories, previously known as Argus Science Inc., Bedford, MA) long-range LRO eye-tracker camera and ET7 software system. Eye-calibration was performed before the resting and task scans, and pupil size was normalized for each individual and scaled to the machine units. Stimulus images were recorded along with the eye image, and simultaneous pupil detection and eye-localization were handled automatically by the ET7 software. Sampling rate was set to 120Hz. To detect eye-blinks, we used ET Analysis (version 184.108.40.206) software’s default settings where blink duration was set to a minimum of 0.1s and a maximum of 0.4s, and to a minimum pupil diameter of 25 units. Temporary pupil disappearances within these values were recorded as “blinks”. Any other interval of data where the pupil was not detected was classified as “signal loss”. We simultaneously registered stimulus event onset times. We tested the reliability of the software performance by eye-balling the pupil size/blink moments and the eye-videos provided by the ASL software on ten randomly chosen blocks from five subjects. For all the subjects we also checked the reliability of the pupil detection algorithm and excluded subjects with major problems as discussed below. The software output accurately reflected the manually recorded blinks.
For all scans, we calculated blink rates estimated as the number of blinks per number of seconds that the eyes were open per each scan run. We standardize this measure across subjects and used it as a covariate in the SPM model.
Task Related Blink Analysis
For the choice-reward task, we calculated the blink densities (blink related eye-closures divided by the number of samples) around the response and feedback stimuli between -3 to +4 seconds. We also calculated the correlations between response times and immediate blink onset latencies after Go stimuli.
Using AFNI 37, timing differences in slice acquisition were corrected, and each volume was aligned to the first volume using rigid body translations and rotations. To remove physiological and motion-related artifacts, a dual-mask spatial independent component (sIC) analysis was applied to the motion and slice-time corrected functional data for each subject 38. The resulting components were classified by a machine learning algorithm based on a simple set of spatial features (including out-of-brain ratio, scattering degree etc.). Brain activity was reconstructed by projecting the selected sICs back to the brain space. The results of the algorithm were monitored and checked by a human expert. The functional and anatomical images of each participant were co-registered and normalized to the stereotaxic space of the Montreal Neurological Institute (MNI) using SPM12. The functional images were resampled to a voxel size of 3 x 3 x 3 mm and smoothed with an isotropic 6-mm FWHM Gaussian kernel. Block-level variance fluctuations were removed in a separate SPM model. Finally, low-frequency (< 1/128 Hz) scanner drifts were filtered with a high-pass filter. The resulting residual activity was then used for cross-correlation analysis with the blink regressors extracted as described below.
Blink related brain activation analysis
For each scan series, blink activity extracted from the eye-tracking analysis was convolved with the canonical Hemodynamic Response Function (HRF) in SPM12 (2-gamma function), and then normalized (z-scored). These time series were used for cross-correlation analysis with voxel-wise brain activity in temporal TR (2.13s) lags from -3 to +3 TRs, and the correlations were then transformed to Fisher's z-values.
Data and participant exclusion criteria
Due to the experiment’s complexity (i.e., setting up eye-tracker in an fMRI environment), we expected some well-known problems to emerge such as i) technical problems in eye-tracking quality (i.e., calibration problems, distortions in pupil detection), and ii) non-compliant participants who did not perform the task properly or closed their eyes too often (particularly during the resting scan). Thus, we excluded any runs/sessions where the task and rest eye data had technical problems, such as long eye-closures (i.e., eye closures >40% of the time) and weak corneal contrasts that interfered with blink detection. This resulted in exclusion of two participants, and one session from each of four participants. Thus, for the final analysis, we had twenty-four sessions collected from fourteen participants with a total of seventy-two task runs and forty-one resting runs available.
Regions of Interest (ROI) Analyses in the Ascending Arousal Network (AAN) and Thalamus
We computed correlation activity plots for regions of interest ROI extracted from nine brainstem nuclei from the AAN (https://www.nmr.mgh.harvard.edu/resources/aan-atlas) using 1mm, MNI152 template 22,39, that comprised the dorsal raphe (DR), median raphe (MR), median raphe frontalis (MRF), locus coeruleus (LC), parabrachial complex (PBC), pontis oralis (PO), periaqueductal grey (PAG), pedunculopontine nucleus (PPN) and ventral tegmental area (VTA). AAN masks are shown in Supplementary Methods. We averaged the z-values for the serotonergic nuclei (MRF, MR, DR) and for the cholinergic nuclei (PO, PPN) and computed separate values for VTA (dopaminergic), LC (noradrenergic), PBC (glutamatergic) and PAG 40. We extracted an additional spherical ROI (3mm diameter) centered at MNI coordinates xyz of 0, -15,-12, which has been used as the location of the VTA ROI by some studies 41. Thus, we wanted to compare findings in VTA when using coordinates from the AAN atlas versus this ROI location, which is slightly anterior and superior to the AAN ROI. For the thalamus we extracted 7 ROIs from the atlas defined by Behrens et al. 28,42 which are found to show exclusive white matter connectivity probabilities with major cortical areas (i.e. motor, premotor, frontal, occipital, temporal, parietal and sensory cortices; masks are shown in Supplementary Methods).
Whole Brain Analyses
For whole brain analyses of blinks associations with brain activity we used SPM12 paired t-test, using the averaged correlation coefficients per subject (N=14) for Task and Rest conditions and included blink rate as covariate (leading to df=12). Normalized (i.e., Fisher’s z-transform) correlations (between BOLD and blink regressor) values of all available resting and task blocks were used in the statistical analysis separately for each lag.
We first conducted a cluster-uncorrected analysis for visualization purposes in brain stem regions (voxel t-threshold t>1.8, k>1, 1mm isotropic, p<.05, liberal approach). In the second analysis we conducted traditional FWER corrected cluster-based analysis (3mm isotropic images) to identify significant clusters across large brain areas. The resolution element size, i.e., ‘resel’, of the images were minimum of 22 voxels leading to around 2985 resels in a group mask, and FWHM was 8.3mm, 8.4mm, and 8.5mm with voxel height thresholding of p<0.002 (corresponding to t=3.54). Clusters (k>30) below pFWER-clus < 0.05 are shown in Figs. 4 and 5 and summarized in Table 1. Detailed cluster maps are given in Supplementary Data 1, 2 and 3 and detailed table of results is given in Supplementary Data 4.
Cerebellar ROI Analyses
Because of the extensive blink associated cerebellar activation signals detected by the whole brain analyses we also extracted 17 ROIs from the cerebellum using the functional connectivity atlas defined by Buckner et al.43, which possess distinct functional connectivity patterns with cortical subregions (See also 44). To simplify our analysis we averaged the values from most of these ROIs to represent their connectivity to 7 main networks defined by the spatial proximity of the target cortical areas (i.e., fronto-parietal attention, DMN, anterior cingulate cortex (ACC), motor, insula, basal forebrain, visual networks). Details of the subregion grouping are given in Supplementary Table 1.