NeuroImage 106 (2015) 414–427

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Coupling between pupil fluctuations and resting-state fMRI uncovers a slow build-up of antagonistic responses in the human cortex Dov Yellin, Aviva Berkovich-Ohana, Rafael Malach ⁎ Department of Neurobiology, Weizmann Institute of Science, 76100 Rehovot, Israel

a r t i c l e

i n f o

Article history: Accepted 16 November 2014 Available online 20 November 2014 Keywords: Pupil dilation Resting-state spontaneous fluctuations Default mode network Locus coeruleus

a b s t r a c t Even in absence of overt tasks, the human cortex manifests rich patterns of spontaneous “resting state” BOLDfMRI fluctuations. However, the link of these spontaneous fluctuations to behavior is presently unclear. Attempts to directly investigate this link invariably lead to disruptions of the resting state. Here we took advantage of the well-established association between pupil diameter and attentional gain to address this issue by examining the correlation between the resting state BOLD and pupil fluctuations. Our results uncover a spontaneously emerging spatiotemporal pupil-BOLD correlation whereby a slow buildup of activity in default mode areas preceded both pupil dilation and wide-spread BOLD suppression in sensorimotor cortex. Control experiments excluded a role for luminance fluctuations or fixation. Comparing the pupil-correlated patterns to activation maps during visual imagery revealed a substantial overlap. Our results indicate a link between behavior, as indexed by pupil diameter, and resting state BOLD fluctuations. These pupil dilations, assumed to be related to attentional gain, were associated with spontaneously emerging antagonism between fundamental cortical networks. © 2014 Elsevier Inc. All rights reserved.

Introduction Although originally considered to be meaningless noise, fluctuations of brain activity during rest, termed “resting state” or spontaneous fluctuations, have recently become a focus of intense research. These ultraslow fluctuations of neuronal activity were first revealed in anesthetized animals (Arieli et al., 1996) and in human BOLD fMRI in the motor system (Biswal et al., 1995). Numerous subsequent studies have validated their significant organization and link to neuronal activity (Greicius et al., 2003; Nir et al., 2008). In particular, these fluctuations have provided new insights into neurocognitive biases and brain pathologies (Agosta et al., 2012; Harmelech and Malach, 2013; Salomon et al., 2012; Wang et al., 2012). However, the precise link of the spontaneous fluctuations to human behavior is still unclear. A number of studies attempted to address this issue by instituting various “free” paradigms. For example, Hesselmann et al. presented subjects with ambiguous visual stimuli at random times and were able to show a small component of anticipatory increase in activity of areas selective to the percept, suggesting a role for spontaneous fluctuations in biasing perceptual decisions (Hesselmann et al., 2008, 2011). Similarly, studying the Libet “free decision” paradigm, Schurger et al. provided an elegant interpretation of the EEG readiness potential as a component of spontaneous fluctuations (Schurger et al., 2012). Finally, recording single unit activity under a “free recall” paradigm, Gelbard-Sagiv et al. discovered a slow

⁎ Corresponding author. E-mail address: rafi[email protected] (R. Malach).

http://dx.doi.org/10.1016/j.neuroimage.2014.11.034 1053-8119/© 2014 Elsevier Inc. All rights reserved.

anticipatory increase in firing rate, again suggesting a contribution of spontaneous fluctuations to free recall (Gelbard-Sagiv et al., 2008). However, a fundamental methodological conundrum in all these attempts to directly examine “spontaneous” behaviors is that they involve an externally specified task. Such a controlled task is likely to interfere with the spontaneous nature of the resting state fluctuations. Even attempts to interrogate individuals at random intervals during rest (Andrews-Hanna et al., 2010) are likely to introduce expectation effects and thus modify subjects' mental state. An alternative approach to studying the potential behavioral links of resting state fluctuations is to measure a behavioral marker, covertly, without subjects' awareness, while they are in the resting state. Indeed, in a previous study of this kind, Ramot et al. were able to demonstrate a significant link between spontaneous BOLD fluctuations and slow spontaneous eye movements that subjects invariably make during rest with their eyes closed. Interestingly, despite the significant link to BOLD fluctuations, subjects were completely unaware of producing such slow eye movements (Ramot et al., 2011). This study demonstrated the feasibility of using covert measures of behavior to uncover a link between behavior (in this case oculomotor drift) and resting state fluctuations. However, a substantial limitation of this study was a lack of any known link between the slow eye movements and any kind of cognitive processing. Here we attempted to advance on this issue by opting to examine the possible link between pupil diameter and spontaneous fluctuations during rest. In contrast to the slow spontaneous eye movements, the neuronal and cognitive processes that control pupil diameter have been extensively studied in humans (Goldwater, 1972) as well as in behaving monkeys (Rajkowski et al., 1994). Importantly, an extensive

D. Yellin et al. / NeuroImage 106 (2015) 414–427

literature has established a clear link between pupil diameter, arousal and attentional gain (Alnaes et al., 2014; Bradley et al., 2008; Bradshaw, 1967; Hess and Polt, 1964; Kahneman, 1973; Naber et al., 2013; Wierda et al., 2012). Furthermore, the neuronal mechanisms that link pupil diameter and attentional gain have been extensively studied in animal models. This research resulted in a theoretical framework (Aston-Jones and Cohen, 2005) termed the ‘Adaptive Gain Theory’ (AGT), pointing to the role of the locus coeruleus (LC)-norepinephrine (NE) system (Sara, 2009) in modulating cortical network processes. Put in simple terms — AGT links LC activity and its concomitant change in pupil diameter with brain-wide modulation of norepinephrine levels. Norepinephrine in turn has been amply documented to be related to cognitive alertness and arousal (Hurley et al., 2004). The theory posits the existence of cortical neuronal accumulators, whose threshold crossing leads to gain modulation (Aston-Jones and Cohen, 2005). Interestingly, while earlier attempts to show a direct link between pupil diameter and LC activity in the human brain using fMRI were unsuccessful (Astafiev et al., 2010)), recent studies using highly engaging tasks, have since revealed such link (Alnaes et al., 2014; Murphy et al., 2014). Intriguingly, during rest, even under conditions of steady illumination, the pupil has been found to undergo slow irregular fluctuations termed “pupil unrest” or “hippus” (Bouma and Baghuis, 1971; Stark et al., 1958). Here, we hypothesized that these apparently random slow fluctuations in pupil diameter during rest are in fact linked to specific network BOLD activations that emerge spontaneously in the resting state. Importantly, given the established connection of pupil diameter to attentional gain, such a link could suggest a role for attentional modulation in the emerging BOLD fluctuations during rest. Furthermore, the specific networks that are positively and negatively linked to peaks of pupil diameter may suggest the content of the spontaneously emerging cognitive processes associated with the pupil dilations during the resting state. We set out to test this prediction by concurrently measuring pupil diameter and BOLD signal when subjects rested with their eyes open. Our results reveal that the spontaneous pupil and BOLD fluctuations are indeed robustly correlated. Aligning the spontaneous BOLD responses with peaks of pupil dilation uncovered a slow buildup of antagonistic responses between the default mode network (DMN) and sensorimotor cortex.

Materials and methods Participants In all, 43 subjects took part in the study (mean age = 31 years, age range = 26–49 years, 19 females). Among them, 22 participated in the main experiment. All subjects had normal or corrected to normal vision, and no history of neurological disorders. Subjects provided written informed consent to a protocol that was approved by the Tel-Aviv Sourasky Medical Center, IRB ethics committee.

Experiments The study consisted of 6 separate experiments, the order of which was counter-balanced between subjects. In all variations of rest experiments reported in this study, subjects were instructed to relax while maintaining their eyes open, during an 8 minute scan. The continuous tracking of participants' eyes allowed verifying their state of wakefulness during the experiment, however participants showing excessive amount of artifacts in pupil data (N 5% of time-course) were excluded. In experiments that included a fixation point, subjects were asked to maintain their gaze on it, and in the absence of a fixation point they were instructed to keep their gaze straight forward. The scanner room was kept dark except for the gray screen (luminance = 90 cd).

415

‘Rest-fixation’ (main) experiment 22 subjects (mean age = 29 years, age range = 26–48 years, 13 females) were instructed to rest while fixating a small (b 10) red dot at the screen's center. The data of two subjects was discarded from subsequent analyses — one due to excessive head movements during the fMRI scan and the second due to excessive pupil data artifacts. ‘Rest-no-fixation’ experiment A subset of 12 subjects from the main (rest with fixation) experiment (mean age = 31 years, age range = 26–48 years, 7 females) received similar instructions as in the main experiment, except for the request to maintain their gaze at the center of a blank gray screen (no fixation point). The data of one subject was discarded from further analyses due to excessive head movements during fMRI recordings. ‘Alternating luminosity’ experiment A subset of 15 subjects from the main experiment (mean age = 28 years, age range = 26–33 years, 10 females) received similar instructions as in the rest-no-fixation experiment. However, without prior notification to subjects, screen luminosity was modulated continuously during the scan, in a 1/f irregular fluctuating manner, replicating changes during an actual rest-fixation pupil diameter time-course (luminance = 65–120 cd). ‘Eyes-open imagery’ experiment A subset of 15 subjects from the main experiment (mean age = 28 years, age range = 26–33 years, 11 females) underwent an 8 min block-design control experiment (20 s of imagery, followed by 10 s of rest). During imagery blocks, subjects were instructed (via short (b1 s) auditory command) to imagine the view of scenes involving exceptionally high or low illumination levels (e.g. viewing a snow covered mountain in a sunny day or entering a dark cinema, respectively), while fixating a small (b1°) red dot on a uniform screen. The data of one subject was discarded from further analyses due to excessive amount of artifacts in pupil data. ‘Eyes-closed' imagery’ experiment 23 subjects (two of which participated in the main experiment; mean age = 33 years, age range = 26–49 years, 8 females) underwent two 10 min block-design imagery experiments. During imagery blocks, subjects were instructed (via short (b 1 s) auditory command) to imagine with closed eyes either navigation along a familiar route (from home to work), or shelves full of familiar tools (either kitchenware or gardening tools). These 21 s blocks alternated with eyes-closed 12 s rest blocks, as well as other eyes-closed language-related conditions, of irrelevance to the current report. Localizer experiment 18 subjects underwent an external “visual localizer” experiment to identify cortical regions that preferentially responded to various categories of visual stimuli. A standard block design was implemented with the following four different categories: blocks of common man-made objects, faces, buildings and geometric textures (image size was approximately 7° ∗ 7° in visual angle). During each block, nine items per category were presented (800 ms per stimulus; 200 ms inter-stimulus intervals). Blocks were separated by 6 s blank-fixation periods, and each condition was repeated over five blocks in pseudo-random order. To maintain attention, subjects were asked to perform a “one-back” memory task (Levy et al., 2004), that is, indicate by a button press when the current stimulus was identical to the one presented just before it. Stimuli Stimuli were generated on a PC (2.4 GHz, 3 GB RAM, NVidia GeForce 8500GT) and projected via an LCD projector (Epson PowerLite 74c) with

416

D. Yellin et al. / NeuroImage 106 (2015) 414–427

refresh rate of 60 Hz onto a tangent screen, viewed by subjects via a tilted mirror. All presentation used in-house scripts implemented using the Cogent toolbox (http://www.vislab.ucl.ac.uk/Cogent/) for MATLAB (Mathworks). Auditory input was delivered via MR compatible head phones (MR Confon GmbH Magdeburg, Germany). Stability of projection The existence of fluctuations in projector illumination was assessed in a dark room by means of a photodiode (PM100, Thorlabs), which was mounted on the central back side of the fMRI facility's projection screen (at standard MRI setup). In order to monitor luminosity fluctuations, the photodiode output was tracked using a standard oscilloscope. Testing was performed under conditions of constant gray screen (luminance = 90 cd) and through gradual shift of grayscale from black to white (luminance = ~10–150 cd). Besides fast oscillations (~200 Hz) consisting of short (b0.5 ms) spikes, presumably due to the projector's backlight flicker, no other oscillations in the level of illumination were noted. MRI setup The scans were performed on a 3 T MRI scanner (Tim Trio, Siemens), equipped with receive only 12-element head matrix coil (Siemens), at the Weizmann Institute of Science, Rehovot, Israel. A whole-brain T1 anatomical image was acquired for each subject to facilitate the incorporation of the functional data into the 3D Talairach space (3D MPRAGE sequence, TR = 2300 ms, TE = 2.98 ms, TI = 900 ms, flip angle = 9° 1 × 1 × 1-mm voxel). Functional images were obtained by T2*-weighted gradient echo planar imaging (EPI) sequence (TR = 2000 ms, TE = 27 ms, flip angle = 75°, FOV 240 mm, matrix size 80 × 80, scanned volume — 32 axial slices of 3 mm thickness (no gap, 3 × 3 × 4-mm voxel), ACPC). Eye-tracking acquisition and preprocessing Pupil diameter was acquired using a noninvasive MR-compatible infrared “Eyelink-1000” (SR Research, Osgoode, ON, Canada) eye tracker, at a sampling rate of 500 Hz. The output of the Eyelink system (providing a scaled index estimate of pupil diameter) was first corrected for periods of missing pupil data (due to blinks or other data acquisition issues). Correction was achieved, first by removing 100 ms before the onset and 100 ms after the offset of each missing data occurrence (due to tendency for artifacts within this range), and then an interpolation of the gap via an inverse-distance weighted method (Howat, University of Washington, 2007). Scans in which the amount of missing pupil data was larger than 10% (typically due to excessive blinking) were excluded from data analysis. The time course was then z-normalized and resampled (using the standard Matlab method) to fMRI TR resolution (0.5 Hz). Implementation of bandpass filters was adapted from EEGLab (Delorme and Makeig, 2004). FMRI data analysis fMRI data were analyzed with the “BrainVoyager” software package (Goebel, 2000) and with complementary Neuroelf toolbox (www. neuroelf.net) code, enclosed within in-house MATLAB (Mathworks) implementation. The cortical surface was normalized to a Talairach coordinate system (Talairach and Tournoux, 1988) and reconstructed for each subject from the 3D-spoiled gradient echo scan. The obtained correlation or activation maps were superimposed on the unfolded cortex. Preprocessing of functional scans included 3D motion correction, slice time correction and filtering out of low frequencies up to two cycles per experiment (slow drift). Subject scans showing sharp head movements (N1 mm) were excluded from the analysis. In order to correlate pupil size with the fMRI BOLD signal, we convolved the preprocessed pupil diameter time-course with a typical hemodynamic response

function (Boynton et al., 1996; Ramot et al., 2011), to create a predictor for standard GLM analysis. The HRF was corrected for the delay in pupil response by subtracting 1 s from the canonical hemodynamic 6 s lag to peak. The analysis was performed independently for the time-course of each individual voxel. After computing the coefficients for the regressor, a student's t-test was performed. In calculating p values, correction for the auto-regression factor was taken into account, due to possible statistical dependencies in the hemodynamic response. Regressors based on white matter and cerebrospinal fluid fMRI signal, as well as on head movement, were included in the model. Multi-subject statistical analysis was based on a random-effect GLM which was FDR corrected whenever the number of analyzed instances was larger than 15, and clustersize corrected via a Monte Carlo simulation (AlphaSim by B. Douglas Ward) otherwise. Regions of interest: definition and analysis Two areas served as the main regions of interest in analysis of this study, based on maximum positive (pIPL) and negative (EVA, approximately centered in V3) pupil-BOLD signal correlation. To avoid circularity, ROIs defined in the main experiment were used only in circumstances — such as band-wise correlational analysis which were unbiased by the method of ROI selection. In Fig. 3, pIPL and EVA ROIs defined in the main experiment were analyzed at different frequency bands, and were cross-correlated with the pupil time-course. These ROIs were also analyzed using cross-spectrum (an estimate of mutual cross power spectral density of two time series, i.e. the distribution of power per unit frequency shared across both time series) and magnitude squared coherence (an estimate of frequency-domain correlations between two time series, indicating how well they correspond at each frequency). When comparing effects related to the strength of pupilBOLD correlations, we avoided circularity by alternating sources of ROI selection between experiments. Thus, in comparison of pupil-BOLD correlations in rest-no-fixation vs. rest-fixation (Fig. S5a) and in eyes-open imagery vs. rest-fixation (Fig. S6a), ROIs defined in one experiment were used for correlation measurement in the complementary experiment and vice-versa, to avoid bias. As for the definition of visual and default-mode networks and ROIs, standard localizers were used. Visual cortex regions of interest (ROIs) were defined individually per each subject via a visual localizer procedure (Levy et al., 2004), as clusters of larger than 200 contiguous voxels (obtained at p b 0.05, FDR corrected). In the resulting map of visual activations (Fig. S1a) yellow–orange regions represent areas that were preferentially activated by visual stimuli. The parahippocampal place, fusiform-face and lateral-occipital areas (PPA, FFA and LO respectively) were defined using the contrasts ‘houses N faces’, ‘face N houses’ and ‘objects N textures’ respectively. Default mode ROIs, including the medial prefrontal cortex, inferior parietal lobule and the junction of the precuneus and posterior cingulated (MPFC, IPL and PCUN respectively), were defined using a method combining in-subject task and rest scans (Salomon et al., 2013). According to this method, subject specific task-negative areas were marked via deactivations during visual stimulation. Then, for each subject, the DMN was defined separately, once using the tasknegative pIPL as a seed, and once the task-negative PCUN as a seed. Thus, average BOLD signal time-course of the pIPL and PCUN seeds during the rest scan was extracted separately and used in a random effect multi-subject analysis of functional connectivity (FC). The resulting significant (r N 0.5, FDR corrected) correlation map was defined as depicting the DMN (Figs. 2b and S3c) and was used to examine the overlap with ROIs defined in the rest-fixation experiment. Finally, a separate analysis was applied for the anterior and posterior parts of the IPL. To that end, an IPL ROI, based on the PCUN seed FC DMN network was sampled using the 200 most anterior and 200 most posterior voxels, defined as aIPL and pIPL areas respectively, in the correlational analysis presented in Figs. S3c–d.

D. Yellin et al. / NeuroImage 106 (2015) 414–427

417

Granger causality analysis

Dynamics of cortical activation associated with pupil dilation movie

We implemented Granger causality (GC) analysis to evaluate linear directed influence between ROIs (which included the conjectured signal from LC, as reflected by pupillary movement). According to the GC model, a time series x Granger-causes y if knowing the past of x helps to predict the future of y better than using the past of y alone. With this in mind, GC analysis, performed using the Granger Causal Connectivity Analysis toolbox (Barnett and Seth, 2014; Seth et al., 2013), was applied to the BOLD time-courses of the pIPL and EVA ROIs, as well as to the corrected and resampled (to 0.5 Hz) pupil time-course. Using methods of this toolbox, we confirmed that our data were covariance stationary (by the augmented KPSS test). Additional preprocessing steps included linear detrending and demeaning. The optimal model order selected using the Bayesian Information Criterion was 6, corresponding to a lag of 12 s. Model validity was confirmed by Durbin–Watson residual whiteness test and the statistical percent consistency test. GC F-statistics were obtained for every ROI pair in both directions and were pooled across subjects.

A movie was constructed from single pupil-BOLD cross correlation map snapshots at different lags (Movie S1). The movie is composed of 4 cycles altogether — two from single subject analysis results shown on a three-dimensional folded map and two from a multi-subject analysis on a cortical unfolded map (N = 20). A timer displays elapsed time relative to the estimated threshold crossing in pIPL, causing LC firing, starting 6 s before and finishing 12 s after each such event.

Shuffled control To create the shuffled control analysis (Fig. S5b), predictors were shuffled randomly between scans of the rest-fixation experiment. A multi-subject map was created by applying a random-effect GLM, in which each subject's dataset was fitted with a design matrix based on another subject's pupil predictor (comprising a single shuffling result). To check the significance of correlation values in key ROIs, a permutation test was conducted. In each instance of the permutation test, all pupil time-course predictors were randomly shuffled between subjects, while the labeling of BOLD time-courses was kept unchanged. The mean pupil-BOLD signal correlation was then calculated per each of the ROIs across subjects. To avoid circularity, ROIs were defined based on restno-fixation experiment. This process was repeated for 10,000 iterations, and the resulting distribution is displayed in the histogram.

Pupil-BOLD cross-correlation approximation A model for approximating the obtained mean pupil-BOLD crosscorrelation from estimated neuronal activity was constructed in stepwise manner. Since the problem is ill-defined (with various degrees of freedom open for arbitrary setting), we began with a highly simplified set of assumptions about the neuronal activations. These neuronal profiles (in pIPL, EVA and LC separately) were then convolved with known hemodynamic and pupil transfer functions (Boynton et al., 1996; Wierda et al., 2012). The predicted pupil and BOLD time courses were then cross correlated as in the experimental results and the two simulated correlation functions were compared to the experimentally derived ones. The simulation was modulated in two additional steps in an attempt to optimize the fit. The fit between the simulation and experimental data was quantified using sum of least squares difference, which was applied both to the mean result as well as for each subject separately. A nonlinear optimization gradient descent algorithm (Matlab ‘lsqcurvefit’ method) was applied over the free parameters in each step, to obtain this minimized distance fit. Significance of fit differences across different simulation steps was assessed using ANOVA and paired t-test (between steps). It should be noted that even better approximation of the actual results was achieved using triangular-like slopes for the neuronal estimated of all three areas at step 3 (not shown). However, this was left out of current scope for the sake of simplicity.

Contribution to overall spontaneous fluctuations Results To perform this analysis (Fig. 7), we defined two ROIs in areas that showed minimal pupil-BOLD correlation during rest-fixation (superior post-central sulcus (sPCS) and dorsolateral prefrontal cortex (DPL)). The size of these ROIs was set to be equivalent with those of the significantly correlated ROIs in pIPL and EVA (i.e. 400–600 voxels) and their mean BOLD activation was also selected to be of similar order (mean(std) in: pIPL — 823(152), EVA — 769(87), sPCS — 881(84) and DPL — 852(69)). Note that due to known B0 and B1 field inhomogeneities in BOLD sensitivity, trying to reach an exact match of mean activation would be impractical. For each ROI, a power-spectrum was computed on the mean time-course within subject. Mean power and its standard error were then established across subjects.

Post-scan interview At the end of each scan, subjects provided a verbal rating response to a set of questions aimed at assessing their experience and ability to perform the task. All questions were rated on a scale of 1–5 (1 being the minimum score and 5 the maximum score). Analysis of variance (ANOVA) for independent samples was used to compare subjects' answers. Post-hoc two-sample t-tests were performed to investigate significant main effects in more detail. For questionnaire content, mean and standard deviation of subjects' replies in relevant experiments, and a statistical analysis pertaining to these replies, see Table S1.

Properties of the rest-state pupil diameter time-course In the main experiment of this study (‘rest-fixation’, see Materials and methods), 22 subjects were scanned while continuously fixating a small dot over a uniform screen for 8 min. Subjects rested without any explicit cognitive task. Following the scans, subjects were interviewed via a short questionnaire (Materials and methods). Subjects' reports indicated high degree of mind wandering during rest-fixation and low attention to the blank screen (Table S1). During the entire scan duration, pupil diameter was measured using an MR compatible eye-tracker (Materials and methods). Fig. 1a depicts a typical trace of pupil diameter changes during the scan. The raw pupil time course, recorded at 500 Hz, was blink corrected, z-normalized and resampled to fMRI acquisition resolution of 0.5 Hz. Power-spectrum and autocorrelation analysis of the raw pupillary responses are depicted in Figs. 1c–d, respectively. The pupil power-spectrum function followed an approximate 1/f slope (of −2.05 +/− 0.12 SEM) up to ~100 Hz. Autocorrelation analysis uncovered no common periodic frequencies. To assess the response dynamics of the pupillary reflex we conducted a separate control experiment in which the screen luminance was changed while the pupil diameter was monitored (see Materials and methods). Cross-correlating luminance and pupillary diameter timecourses revealed the previously described response latency of ~ 1 s (both to luminance change and attentional events) (Ellis, 1981; Wierda et al., 2012), as shown in Fig. 1e.

418

D. Yellin et al. / NeuroImage 106 (2015) 414–427

Fig. 1. Rest-fixation experiment — pupil diameter time-course analysis. Example from a single subject's pupil time-course during the main experiment — (a) raw pupil diameter signal in gray and blink-corrected is shown in blue. (b) Same subject's HRF-convolved pupil diameter time-course shown in red. (c, d) Power spectrum and auto-correlation analysis of pupil diameter changes, respectively (N = 20) — average shown in blue overlaid on individual subjects data in gray. (e) Cross-correlation of pupil and luminosity time-courses taken from the alternating luminosity experiment.

Pupil-BOLD time-course correlations during rest with fixation To examine whether the spontaneous fluctuations in pupil diameter were correlated with BOLD fluctuations, we convolved the resampled pupil time-course with a standard hemodynamic response function (Boynton et al., 1996) and used this as a predictor (see Fig. 1b) in fMRI GLM analysis (Ramot et al., 2011). To compensate for the abovementioned delay in pupil response, the canonical ~6 s to peak hemodynamic response lag was shifted back by 1 s. Results revealed a highly consistent map of positive and negative correlations to the pupildiameter predictor throughout the cortex. Fig. 2 depicts the pattern of this pupil-BOLD signal correlation for a single representative subject (Fig. 2a) and for the entire 20 artifact-free subject group (Fig. 2b). Regions of significant pupil to BOLD negative-correlations were widespread in regions of the sensorimotor ‘extrinsic’ system (Golland et al., 2007; Ramot et al., 2011). Specifically, robust anti-correlations were found in visual, motor and somatosensory regions. Several positively correlated areas were noted as well, anatomically identified as posterior inferior parietal lobule (pIPL, appearing most consistently across subjects), the junction of posterior cingulate cortex-precuneus (PCUN) and medial prefrontal cortex (MPFC). To identify the relationship of this pattern to known visual areas as well as to intrinsically-oriented DMN, we implemented visual and DMN functional localizers (Materials and methods and Fig. S1). As can be seen in Fig. 2b, areas activated by visual stimuli (in purple) included

the area of highest negative pupil-BOLD correlation, but not the positively correlated areas (pIPL, PCUN and MPFC) which fell within the boundaries of the DMN (in cyan). For ROI analysis, two regions exhibiting the strongest and most consistent positive (pIPL) and negative (early visual, EVA) pupil-BOLD correlation (Materials and methods and Fig. S2) were selected. Using these ROIs we examined to what extent the pupil-BOLD correlation was evident in specific frequency bands. To that end we calculated the cross-spectrum and magnitude squared coherence relationships (see Materials and methods) between the pupil trace and average timecourse of each ROI. The results across 20 subjects of the main experiment, shown in Figs. 3a–b, demonstrate that for both negative and positive correlations, slow frequencies (b 0.05 Hz) dominated the common spectrum's power and mutual correlation. Figs. S3a–b depict rest-fixation pupil-BOLD correlation maps based on high (N 0.2 Hz) vs. low (b0.05 Hz) bandpass filters (Materials and methods), correspondingly. The maps clearly reveal dominance of lower frequencies in the resulting pattern of correlations. An important issue concerns the component of spontaneous BOLD fluctuations actually explained by the pupil diameter changes. Examining this issue in the coherence analysis, revealed that at cortical regions manifesting the highest positive and negative correlations, the coherence level at the optimal frequencies (b0.05 Hz) reached values of: 0.36 +/− 0.042 and 0.48 +/− 0.059 (mean & standard error), respectively, explaining up to ~25% of the variance.

D. Yellin et al. / NeuroImage 106 (2015) 414–427

419

Fig. 2. Spatial topography of pupil-BOLD signal correlations in rest-fixation experiment. (a) Single subject correlation map, projected on an inflated (top) and unfolded cortex (below). (b) Multi-subject pupil-BOLD correlations on an unfolded map. Random effects analysis, corrected, on 20 subjects. Note that in both maps, sensorimotor areas are negatively correlated and default mode areas are positively correlated to the pupil diameter predictor. Color scale indicates statistical significance. Yellow–orange regions represent areas for which the BOLD signal was positively correlated to the pupil predictor, whereas blue–green regions indicate negative correlations. Purple and green contours delineate the borders of the visual and default networks, respectively CS, central sulcus; EVA, early visual cortex; IPL, inferior parietal lobule; IPS, intraparietal sulcus; LS, lateral sulcus; MPFC, medial prefrontal cortex; PCUN, precuneus; LH, left hemisphere; RH, right hemisphere; A, anterior; P, posterior.

Temporal properties of pupil-BOLD signal correlation So far our analysis focused on the “static” correlation pattern, i.e. the pupil-BOLD signal coupling, taking into account their respective expected delayed response. An interesting question concerns the dynamic evolvement of the correlation between these signals, i.e. how do they develop over time. To address this issue we cross correlated the restfixation pupil diameter and BOLD signal time courses. The results, as depicted in Fig. 3c., indicate that maximal negative correlation in EVA ROI BOLD signal peaked at 5.31 +/− 0.35 s SEM lag relative to the pupil diameter signal, falling within the approximate lag of coactivation with the LC (assumed to drive the pupil dilation (Ellis, 1981)). Interestingly, the maximal positive correlation appeared in pIPL ROI at 2.47 s +/− 0.61 SEM lag between the pupil and BOLD time courses, i.e. preceding the predicted LC activity by ~2.5 s. This difference in peak cross-correlation was found to be significant across subjects (p b 10−4, one-sided t-test, see inset). Other DMN nodes (in PCUN and MPFC ROIs) exhibited early peak cross-correlation as well, although weaker and somewhat deferred relative to pIPL. To examine whether the effect may have been dependent on the manner by which the pIPL

ROI was defined, we compared the effect also with a pIPL ROI defined through FC to the PCUN seed — which may provide a more reliable definition of DMN (see Materials and methods). Comparing the pIPL ROIs defined based on the PCUN seed vs. the task-negative responses in the pIPL revealed no significant difference (paired t-test: p = 0.16). A subdivision of the IPL (based on PCUN seed FC analysis — see Materials and methods) into anterior and posterior ROIs (aIPL and pIPL respectively) revealed a significant difference in the strength of pupil correlation in the latter (paired t-test: p b 0.05). Pupil-BOLD cross-correlation across the entire cortex on an unfolded map, in movie format (see Movie S1) showed a similar pattern of surprisingly early positive and late negative correlations. To better assess the relationship between events associated with these pupil-BOLD correlations, and in view of the unexpected long lags, we applied Granger causality analysis (Seth et al., 2013) (see Materials and methods). Group results shown in Fig. 3d demonstrate significant directionality effects from pIPL onto the pupil (assumed to reflect LC activity, FpIPL → pupil − Fpupil → pIPL = 6.78, p b 0.01) and from the pupil onto EVA (Fpupil → EVA − FEVA → pupil = 13.94, p b 0.01). Interestingly, the interrelation from pIPL to EVA showed a weak trend only (p = 0.13). These

420

D. Yellin et al. / NeuroImage 106 (2015) 414–427

Fig. 3. Pupil-BOLD correlation during rest-fixation — static and dynamic ROI analysis. Pupil-BOLD (a) cross-spectrum and (b) magnitude squared coherence analyses (N = 20) in the posterior inferior parietal lobule (pIPL, in red) and early visual area (EVA, in blue) ROIs. Note the significant contribution of slow frequencies to the common correlation. (c) Multi-subject (N = 20) cross-correlation plots of pupil-BOLD time-courses, for the same subjects, the two ROIs as in panel a (pIPL and EVA) and two additional DMN ROIs — PCUN (in magenta) and MPFC (in green). Inset panel on the left shows the group average of lag to max-correlation in the two main pIPL and EVA ROIs. (d) Granger causality analysis indicating directionality relations between pupil and BOLD time-courses. The histogram depicts the mean connectivity f-statistics across subjects, between all possible edge arrangements of the pupil and the two ROIs (i.e. pIPL to pupil, EVA to pupil, pIPL to EVA and vice-versa directions). The inset figure demonstrates the resulting significant directionality effect, a strong connectivity from pIPL onto the pupil, and from the pupil onto EVA.

findings were further validated via a permutation test (Fig. S4) indicating that the two significant links fell significantly beyond chance level (p N 0.99), as opposed to the other configurations. Thus, both Granger and cross-correlation analyses appear to indicate a cascade of events beginning with pIPL activity and culminating in EVA inhibition (inset). Robustness of the rest-fixation pupil-BOLD signal correlations Could the need to maintain fixation have led to the observed correlation patterns? To examine this issue we conducted a control experiment in which subjects were asked to view a blank field with no fixation (rest-no-fixation). Our results (Fig. S5a) show similar anatomical distribution of positive and negative correlations, albeit somewhat weaker (0.14 +/− 0.021 (SEM) in pIPL and − 0.23 +/− 0.013 in EVA) in no-fixation compared to (.19 +/− 0.012 and 0.27 +/− 0.030, respectively) with fixation. This difference was statistically significant in pIPL (p b 0.05, paired t-test), but not in EVA. It should be noted that subjects' post-scan questionnaire rankings indicated a significant difference between these two conditions as well (Table S1 — question # 2). Another potential confound may have been slow fluctuations in screen luminance — either due to some instability in our projector

device (but see Materials and methods) or because of ocular aspects, e.g. changes in luminance due to the pupillary fluctuations themselves. Furthermore, an efferent copy or feedback from the pupil motor control network itself may have initiated the observed BOLD fluctuations. To explore these potential factors, we examined BOLD responses when screen luminance was artificially modulated in a similar fashion to typical pupil diameter fluctuations (Materials and methods). As depicted in Fig. 1E, the luminance changes were effective in driving pupil responses. The impact of these changes on the BOLD response is shown in Fig. 4a. Evidently, change in luminosity failed to induce significant widespread responses, with the exception of a small region of increased luminance activation (corresponding to narrowing of pupil diameter) in foveal V1–V2 (arrow), in agreement with earlier studies (Rossi and Paradiso, 1999; van de Ven et al., 2012). The map depicted in Fig. 4b shows weaker pupil-BOLD correlations during the alternating luminosity experiment, relative to the main experiment. Yet, regressing out the luminosity change profile restored stronger correlation patterns as apparent in the main experiment (Fig. 4c). To investigate the possibility that the correlation maps were a product of some spurious auto-correlation structure of the pupil and BOLD fluctuations, a shuffle control was conducted. A map constructed by

D. Yellin et al. / NeuroImage 106 (2015) 414–427

421

Fig. 4. Alternating luminosity control experiment — multi–subject analysis. (a) Screen Luminosity-BOLD correlation map. Screen luminosity was modulated similar to pupil diameter changes and used as a predictor of the fMRI BOLD signal (see Materials and methods), shown on unfolded cortical map. Note the minimal impact of the physical changes in luminosity or luminosity-driven pupil changes on the BOLD response, which was restricted to early visual areas. (b) Pupil-BOLD correlation map during rest with alternating screen luminosity. (c) Pupil-BOLD correlation map based on GLM with luminosity set as a confound. In all maps N = 15, random effects, p b 0.05, corrected. Abbreviations — same as Fig. 2.

randomly shuffling subjects' pupil predictors and a complementary permutation test (Materials and methods and Fig. S5B) failed to reveal significant correlations, attesting to the robustness of our results. Regressing out of white matter and cerebrospinal fluid fMRI signal (assuming they may reflect some confounding component of autonomic nervous system activity) had only a weak impact on the obtained correlation maps. Comparison of rest with visual imagery Could the consistent pattern of rest-fixation pupil-BOLD correlations be indicative of some underlying self-initiated cognitive process? To examine this potential relationship we compared them to activation patterns generated during an internally-elicited process, namely a visual imagery task, reported to also give rise to pupil dilation (Reeves and Segal, 1973; Simpson and Paivio, 1966). A block-design ‘eyes-open imagery’ experiment was implemented, in which subjects fixated a uniform field in a similar fashion to their rest-fixation experiment. We asked subjects to imagine episodes in which they were exposed to either intensely bright or very dark surroundings (Materials and methods). The results of this experiment are depicted in Fig. 5, in which we compare patterns of pupil-BOLD correlations emerging during the rest-fixation (Fig. 5a) and the eyes-open imagery experiments (Fig. 5b) with imagery evoked activation (Fig. 5c). Examination of the maps reveals a clear, albeit partial, overlap between the imagery-

related BOLD responses and the pupil-BOLD correlations during both rest and imagery. Intriguingly, a particularly robust common effect was the significant negative correlation of early visual areas — evident both in the rest-fixation condition (blue contours) and during imagery. As for positive correlations (red contours), they appeared to overlap mainly in the pIPL, while a substantial segregation was apparent in the medial-parietal DMN structures (PCUN and the MPFC). Yet, a major discrepancy between the imagery and rest-fixation maps was evident in high order visual cortex – particularly the LH anterior PPA and POS – which were activated during imagery but showed no significant pupil-BOLD correlations during rest (Figs. 5b–c, white arrows). A direct statistical comparison between the ROI-specific pupil correlations in rest-fixation vs. eyes-open imagery (Figs. S6a–b) showed a trend for stronger correlations in both selected DMN ROIs (pIPL and PCUN) during rest (paired t-test: p = 0.14 and p = 0.17 respectively). Interestingly, the negative correlations in a high-order visual area (PPA) showed a weak trend for stronger activation during imagery (p = 0.24), whereas an opposite trend was observed (i.e. weaker activation in imagery compared to rest) in low-order visual areas (EVA, p = 0.09). Regions exhibiting imagery-related activation and inactivation in the contrast map were generally in close correspondence with the pupil correlated regions. Nevertheless, hemispheric asymmetry was noted in both imagery related maps (showing a stronger effect in the left hemisphere), whereas a more laterally symmetric map was obtained during restfixation (see gray arrow pointing to RH pIPL). No significant BOLD signal

422

D. Yellin et al. / NeuroImage 106 (2015) 414–427

Fig. 5. Visual imagery with eyes open experiment — multi-subject analysis. (a) Map of pupil-BOLD correlation in the main experiment (rest-fixation) for comparison. Contours from positively (red) and negatively (blue) pupil-BOLD correlated areas are highlighted for later comparison with the imagery experiment results. (b) Map of pupil-BOLD correlation during the eyes open imagery experiment. (c) Activation map for the contrast of imagery N rest. In all maps N = 14, random effects, p b 0.05, corrected. Abbreviations — same as Fig. 2. Note the robust overlap in many sensorimotor and DMN areas, and also the non-overlapping regions in rest vs. imagery (indicated by white and gray arrow-heads), such as the right parietal and PPA regions.

difference was observed between the imagery of dark vs. illuminated scenes. Could the eyes-open imagery related responses have been due to the explicit instructions the subjects received, or to the fact that their eyes were open? In order to examine these issues we compared the eyes-open imagery results with ‘eyes-closed imagery’ experiments conducted on a separate group of subjects, in which a similar block design was employed. Due to eye closure, we were unable to track the pupil diameter; however, the eyes-closed control allowed us to compare the BOLD activation with that of eyes-open imagery in two equivalent but slightly modified imagery tasks (navigation and tools) of similar nature (Materials and methods). In Fig. 6, the activation map for eyes-open imagery (Fig. 6a) is compared with corresponding maps of these two eyesclosed imagery control experiments (Figs. 6b–c, respectively). It should be noted that the activity baseline in rest blocks could be different in eyes-open vs. eyes-closed, yet as can be seen, the overlap of BOLD activation in the three distinct imagery tasks was substantial. Nevertheless, differences associated with imagery task specificity can be noted as well; for example, the positive activation in high-order lateral occipital (LO) object-selective area (Malach et al., 1995) that was evident only for the tool imagery task. Furthermore, the level of inactivation in early visual areas was markedly reduced in the eyes-closed condition.

A comparison of the average BOLD signal during eyes-open vs. eyesclosed imagery (see Figs. S6c–d) indicated a trend toward higher and longer activity levels in pIPL and EVA during eyes-closed. However, this could possibly be due to a higher baseline of rest activity during the eyes-open state.

Spontaneous fluctuations in the absence of pupil-BOLD correlations A related issue is whether, during rest-fixation, the process underlying pupil fluctuations was the principal mechanism driving the spontaneous BOLD fluctuations. If this was indeed the case, we would expect the amplitude (power) of BOLD fluctuations in regions that were uncorrelated to the pupil fluctuations to be significantly reduced. To examine this issue, we compared the power of the BOLD fluctuations in cortical regions showing high pupil-BOLD correlations (pIPL and EVA), with others showing minimal correlations (sPCS and DPL, see Materials and methods). There was a small, insignificant, trend toward lower power values in the pupil-uncorrelated areas (Fig. 7), indicating that the BOLD fluctuations during rest-fixation were not confined to the pupil correlated areas, but continued independently throughout other parts of the cortex as well.

D. Yellin et al. / NeuroImage 106 (2015) 414–427

423

Fig. 6. Visual imagery — eyes closed vs. eyes open. Activation maps comparing the results of eyes-open and eyes closed imagery experiments. Contrast of imagery N rest is compared between (a) imagery of variably illuminated scenes with eyes-open (N = 14), (b) imagery of navigation and (c) tools with eyes-closed (N = 23) on unfolded cortical maps. Contours of activated (brown) and deactivated (green) areas based the eyes-open imagery (a) are overlaid on maps of both eyes-closed imagery experiments. All maps use random-effect, p b 0.05, corrected. Abbreviations — same as Fig. 2. Note substantial overlap between the maps, as well as distinct discrepancies (arrows).

Discussion Our study examined the hypothesis that the spontaneous (resting state) BOLD fluctuations can be, at least partially, related to spontaneous modulations in pupil diameter. Examining the relationship between spontaneous pupil and fMRI BOLD fluctuations uncovered two consistent findings. First, there was a robust and highly significant pupilBOLD correlation in specific cortical networks. Second, this correlation pattern showed consistent spatiotemporal dynamics. Interestingly, our findings uncovered a slow cascade of antagonistic response, whereby a slow buildup of positive BOLD activation in pIPL (as well as in additional hubs of the DMN/intrinsic system) anticipated the pupil dilation as well as negative BOLD responses in widespread areas of the sensorimotor, extrinsic system (Chai et al., 2012; Golland et al., 2007; Keller et al., 2013; Weissenbacher et al., 2009). The effect was most pronounced in the pIPL and was not dependent on the manner by which the pIPL ROI was defined — either directly from its task-negative responses, or indirectly through its functional connectivity to the task negative PCUN seed. Interestingly, we observed a significant difference in the correlation between DMN subdivisions. Thus, analyzing separately the anterior and posterior parts of the DMN revealed a significantly higher correlation level in the posterior IPL compared to the anterior part (Fig. S3d). This result is compatible with our recent observation of subdivision using task-design (Salomon et al., 2013). In view of the extensive literature linking pupil diameter to attentional gain and arousal (Alnaes et al., 2014; Bradley et al., 2008;

Bradshaw, 1967; Hess and Polt, 1964; Kahneman, 1973; Naber et al., 2013; Wierda et al., 2012), our results strongly suggest that this spontaneously emerging antagonistic cascade is linked to changes in attentional state and cognitive demand. Our results thus contradict the hippus model of random pupil fluctuations (Bouma and Baghuis, 1971; Stark et al., 1958). They point to a new link between resting state fluctuations and behavior. The pupil-BOLD correlations were dominated by very slow frequencies (b~0.05 Hz), arguing against a contribution of fast behavioral events (such as blinks) to this result. It should also be noted that although the hemodynamic response function is slow, it could not account for the ultra-slow frequencies dominating the correlations, further supporting their link to the ultra-slow spontaneous fluctuations (Nir et al., 2008). The neuronal origins of the pupil-BOLD dynamics What could be the temporal sequence of neuronal events that underlies the observed dynamics of the pupil-BOLD correlations? Results of the luminance-control experiment ruled out subtle fluctuations in screen luminosity during the rest-fixation experiment as the driving force. Furthermore, they argue against the possibility that mere changes in pupil diameter could have caused the BOLD responses e.g. via efferent copy mechanism. In an attempt to get a clearer view of the neuronal processes underlying the dynamics of the observed pupil-BOLD correlations, we took advantage of the known transfer functions linking neuronal firing to

424

D. Yellin et al. / NeuroImage 106 (2015) 414–427

Fig. 7. Analysis of BOLD fluctuations power measured at cortical regions showing high and low coupling to pupil. (a) BOLD-pupil correlation during rest-fixation indicating locations of sampled ROIs. White arrows indicate the location of pupil correlated (EVA and pIPL) and uncorrelated (sPCS and DPL) areas, which were used in current analysis. (b) Histograms depicting the BOLD power in each of the ROIs plotted separately for high and low amplitude pupil fluctuations (large and small respectively). The amplitude of BOLD fluctuations in regions uncorrelated to pupil fluctuations showed a small but insignificant, trend toward lower values.

BOLD (Boynton et al., 1996; Mukamel et al., 2005) on the one hand, and LC activity onto pupil diameter (Wierda et al., 2012; Zylberberg et al., 2012) on the other. We used these functions to approximate, in a highly simplified manner, the experimentally derived pupil-BOLD crosscorrelations (Materials and methods and Table S2). Fig. 8 depicts the results of this step-wise optimization process. As can be seen, assuming the simplest option of a concurrent activation in LC, pIPL and EVA failed to capture the observed dynamics (Fig. 8a). A necessary modification that greatly improved the fit of the simulation was the introduction of a lag between pIPL and LC neuronal activations (Fig. 8b). The obtained lag (of few seconds) was clearly far longer than expected by simple synaptic delays, which are typically on the order of fractions of a second (Fisch et al., 2009; Quiroga et al., 2005). Nevertheless, long duration build-ups have been a well-established aspect of accumulator type mechanisms, extensively studied in the domain of decision making (Cain et al., 2013). Interestingly, long duration anticipatory activity has been reported also in spontaneous types of decision processes, such as volitional perceptual decisions (Hesselmann et al., 2011), free recall (Gelbard-Sagiv et al., 2008) and “readiness” activity prior to spontaneous motor responses (Soon et al., 2008). It has been conjectured, that such long delays may point to the operation of accumulator dynamics also in the case of resting state activity (Schurger et al., 2012). In line with this literature, we have attempted to replace the step-function of

the pIPL neuronal activity in our model with a slow, accumulator-like, build-up of activity — and indeed this led to a closer fit to the experimental results (see Fig. 8c). Thus, a plausible scenario that may explain the dynamics of spontaneous events in our experiment posits a slowly developing spontaneous wave of activation in pIPL (and additional DMN nodes) which eventually crosses a decision threshold (Fisch et al., 2009). This threshold crossing triggers a burst of activity in the LC-NE system (Aston-Jones and Cohen, 2005), leading to an increase in attentional gain, pupil dilation and the observed EVA inactivation. However, it should be noted that the accumulator dynamics suggested by our findings appear to be slower compared to those mentioned in former studies. Whether this is due to the sluggish nature of the BOLD signal, or to a true neuronal dynamics remains to be evaluated. Additional subject by subject analysis of the pIPL signal derivative, preceding significant pupil dilation events, might better reveal the details of these dynamics. However, such detailed analysis may necessitate a better temporal resolution, e.g. derived from intracranial ECoG recordings. The exact anatomical site and pathways, that could carry an accumulator threshold crossing signal from DMN nodes (such as pIPL) to the LC, are yet unknown. Being situated at the junction point of the frontoparietal attention (Wang, 2008) and the default mode networks, the posterior parietal cortex has been suggested to play a role in processes involving both

D. Yellin et al. / NeuroImage 106 (2015) 414–427

425

Fig. 8. Model simulation of the Pupil-BOLD dynamics during rest-fixation. The simulation was built as a series of successive approximations starting from the simplest possible assumptions about the neuronal activity in pIPL, LC and EVA (see Supplementary data). In each step, the neuronal estimates are convolved with established transfer functions to yield the expected response profiles which are then cross-correlated as in the actual rest-fixation experiment. The neuronal profile variables (duration, lag and slope) were optimized at each step, using least-squares minimization to reach a best possible fit under the simplifying assumptions set by the given variables. We start with the assumption of concurrent equal duration neuronal activity in all three brain regions. As can be seen, results converge to highly optimized fit by introducing a lag of few seconds between pIPL and LC (step 2) and introducing gradualaccumulator like- buildup of activity in pIPL (step 3). Statistical analysis across subjects (inset) showed that each additional step led to a significantly better fit (ANOVA yielded F(2, 57) = 21.02 and paired t-test between steps indicated p b 0.001).

cognition and attention, such as time-keeping (Buhusi and Meck, 2005). Furthermore, its involvement in self-initiated thought during mind wandering, which has been shown to also be reflected in pupil size (Grandchamp et al., 2014; Smallwood et al., 2012), should be noted. Thus, it would be plausible to consider that the signal accumulation of pIPL activity may underlie the initiation of mind wandering related events during the resting state. Unfortunately, the current experimental design was not aimed at a deep analysis of this interesting possibility. Implications for the adaptive gain theory Having employed a weakly engaging task paradigm impaired our ability to directly correlate pupil diameter with LC activity via the BOLD signal (Astafiev et al., 2010). Yet, in view of previous demonstrations of a link between BOLD fMRI and LC activity (Alnaes et al., 2014; Murphy et al., 2014) our findings appear to also be compatible with the notion that activity in the LC-NE system underlies the observed pupil-BOLD correlation. Our results suggest that a neural accumulator may play a role in task-free activation of LC (i.e. tonic mode). Earlier work on AGT, in which an underlying neural accumulator was hypothesized (Aston-Jones and Cohen, 2005), suggested its role might reflect a decision tradeoff between exploitation and exploration modes of attention, i.e. whether to continue pursuing a current task or divert to an alternative of higher utility. However, since our experiment was not aimed to address the exploitation-exploration conjecture, this will need to be addressed in future studies. Actual experiments thus far

mostly concentrated on task-positive implications of the theory (testing the phasic mode response of LC), demonstrating that results are consistent with it. Nevertheless, it has already been shown that spontaneous thoughts arise most frequently during high tonic levels of LC activity (Smallwood et al., 2012), a finding well compatible with our findings. In agreement with earlier studies in non-human primates (Kobayashi et al., 2000; Logothetis, 2010) and more recently in human LC using BOLD-fMRI (Alnaes et al., 2014; Murphy et al., 2014), the current findings point toward cortical suppressive effects in sensorimotor cortex being linked to gain modulation during periods of increased tonic LC activity, also marked by spontaneous pupil dilation. Comparing the spatial layout of correlation during rest and imagery While the threshold crossing model can account for the slow cascade observed in our data, it leaves open the spatial aspects, i.e. accounting for the localization of positive and negative pupil-BOLD correlations to DMN and sensorimotor structures, respectively. Here we can turn to the imagery results as a source of potential explanation. Our results confirm previous reports showing a positive correlation between pupil dilation and visual imagery (Reeves and Segal, 1973; Simpson and Paivio, 1966). Furthermore, the pattern of pupil-BOLD correlation as well as the functional contrast imagery N rest during the eyes-open imagery experiment revealed a substantial, albeit incomplete, overlap with that observed in rest-fixation (Fig. 5). The observed pattern of positive and negative activations was not a result of the fact

426

D. Yellin et al. / NeuroImage 106 (2015) 414–427

that the eyes were open during the imagery task — since a similar pattern appeared during the more conventional imagery procedures conducted in darkness with eyes closed (Fig. 6). The DMN has been consistently implicated in self-related and episodic memory functions (Buckner and Carroll, 2007). So, obtaining activity in pIPL and other DMN structures during both rest-fixation and imagery may suggest that similar, internally oriented, cognitive processes were generated during these experiments. The intriguing inactivation of EVA during imagery suggests that an important component of visual imagery is paradoxically related to inhibiting spontaneous activation in early visual cortex. This effect is compatible with prior studies (Daselaar et al., 2010) and with our previous suggestion (Amedi et al., 2005) that visual imagery is complemented by suppression of potential sources of distracting information derived from competing sensory representations. This interpretation is further supported by the tendency of the EVA suppression to be less pronounced during eyes-closed imagery (Fig. 6). Our results are also compatible with the interpretation (Einhauser et al., 2008) suggesting that upon threshold crossing, NE gain modulation enhances dominant representations, while simultaneously reducing the activity in competing ones, including their levels of spontaneous neuronal activity. On the other hand, our results argue against models proposing that imagery is dependent on concurrent positive reactivation of early visual areas (Kosslyn et al., 2001). This could perhaps be explained by the imagery paradigm we employed, relying on the need to initiate lengthy self-related reflection vs. the shorter, snapshot-like, imagery tasks used in other paradigms (Ishai et al., 2000). The possibility of subjects' inability to maintain the state of imagery throughout the block can also not be ruled out. Contribution to overall spontaneous fluctuations While our results show a robust and consistent relationship between BOLD and pupil fluctuations, an important question is to what extent can the spontaneous fluctuations be fully explained by pupil diameter changes. This question relates to two separate aspects. First, within the regions that showed such correlations, how much of the BOLD fluctuations correlated to the pupil diameter changes? The second aspect relates to areas that were uncorrelated to the pupil-diameter fluctuations. Did these uncorrelated areas continue to manifest spontaneous fluctuations? Our results show that a substantial component of the BOLD fluctuations could indeed be explained by the pupil diameter fluctuations. However, it is clear that additional modulatory processes were affecting the spontaneous BOLD fluctuations in these regions beside those directly related to pupil diameter changes. Another important observation in this regard was the lack of significant reduction in the amplitude of BOLD fluctuations in areas located outside the pupil correlated regions (Fig. 7). Our analysis revealed the existence of areas exhibiting weak pupil correlation albeit their spontaneous BOLD fluctuations. These results thus support the notion that BOLD fluctuations are indeed incessant and occur throughout the cortex. Nevertheless, they leave open the question why a subset of these fluctuations (e.g. in pIPL) was coupled to pupil dilation while others were not. If our accumulator hypothesis is correct, these results point to other factors, such as increased baseline shifts in the pIPL rather than increased amplitude of its spontaneous fluctuations as the mechanism leading to the selective threshold crossing in this region. Further research will be needed to explore the validity of these conjectures. Conclusion To summarize, our findings reveal a spontaneously emerging antagonism between fundamental networks of the human cortex that is correlated to pupil dilatations during rest. Thus, our results reveal a link between behavior – indexed by pupil diameter fluctuations – and spontaneous network activations and deactivations during the resting state.

Acknowledgment The authors wish to thank Michal Harel for the help in brain reconstructions, and Nahum Stern, Fanny Atar and Dr. Edna Haran-Furman for their technical MRI facility assistance. We also thank Amos Ariely, Dov Sagi, Alumit Ishai, Oren Shriki and the Malach lab members in general for the constructive discussions; in particular, Tal Golan, Meiri Meshulam and Michal Ramot for reviewing the manuscript and providing many helpful comments. This work was funded by the EU FP7 VERE, EU-Flagship HBP, ICORE program (ISF 51/11) and the Helen and Martin Kimmel award to R.M.

References Agosta, F., Pievani, M., Geroldi, C., Copetti, M., Frisoni, G.B., Filippi, M., 2012. Resting state fMRI in Alzheimer's disease: beyond the default mode network. Neurobiol. Aging 33, 1564–1578. Alnaes, D., Sneve, M.H., Espeseth, T., Endestad, T., van de Pavert, S.H., Laeng, B., 2014. Pupil size signals mental effort deployed during multiple object tracking and predicts brain activity in the dorsal attention network and the locus coeruleus. J. Vis. 14. Amedi, A., Malach, R., Pascual-Leone, A., 2005. Negative BOLD differentiates visual imagery and perception. Neuron 48, 859–872. Andrews-Hanna, J.R., Reidler, J.S., Huang, C., Buckner, R.L., 2010. Evidence for the default network's role in spontaneous cognition. J. Neurophysiol. 104, 322–335. Arieli, A., Sterkin, A., Grinvald, A., Aertsen, A., 1996. Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses. Science 273, 1868–1871. Astafiev, S.V., Snyder, A.Z., Shulman, G.L., Corbetta, M., 2010. Comment on “Modafinil shifts human locus coeruleus to low-tonic, high-phasic activity during functional MRI” and “Homeostatic sleep pressure and responses to sustained attention in the suprachiasmatic area”. Science 328, 309 (author reply 309). Aston-Jones, G., Cohen, J.D., 2005. Adaptive gain and the role of the locus coeruleus– norepinephrine system in optimal performance. J. Comp. Neurol. 493, 99–110. Barnett, L., Seth, A.K., 2014. The MVGC multivariate Granger causality toolbox: a new approach to Granger-causal inference. J. Neurosci. Methods 223, 50–68. Biswal, B., Yetkin, F.Z., Haughton, V.M., Hyde, J.S., 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. Bouma, H., Baghuis, L.C., 1971. Hippus of the pupil: periods of slow oscillations of unknown origin. Vis. Res. 11, 1345–1351. Boynton, G.M., Engel, S.A., Glover, G.H., Heeger, D.J., 1996. Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neurosci. 16, 4207–4221. Bradley, M.M., Miccoli, L., Escrig, M.A., Lang, P.J., 2008. The pupil as a measure of emotional arousal and autonomic activation. Psychophysiology 45, 602–607. Bradshaw, J., 1967. Pupil size as a measure of arousal during information processing. Nature 216, 515–516. Buckner, R.L., Carroll, D.C., 2007. Self-projection and the brain. Trends Cogn. Sci. 11, 49–57. Buhusi, C.V., Meck, W.H., 2005. What makes us tick? Functional and neural mechanisms of interval timing. Nat. Rev. Neurosci. 6, 755–765. Cain, N., Barreiro, A.K., Shadlen, M., Shea-Brown, E., 2013. Neural integrators for decision making: a favorable tradeoff between robustness and sensitivity. J. Neurophysiol. 109, 2542–2559. Chai, X.J., Castanon, A.N., Ongur, D., Whitfield-Gabrieli, S., 2012. Anticorrelations in resting state networks without global signal regression. NeuroImage 59, 1420–1428. Daselaar, S.M., Porat, Y., Huijbers, W., Pennartz, C.M., 2010. Modality-specific and modality-independent components of the human imagery system. NeuroImage 52, 677–685. Delorme, A., Makeig, S., 2004. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. Einhauser, W., Stout, J., Koch, C., Carter, O., 2008. Pupil dilation reflects perceptual selection and predicts subsequent stability in perceptual rivalry. Proc. Natl. Acad. Sci. U. S. A. 105, 1704–1709. Ellis, C.J., 1981. The pupillary light reflex in normal subjects. Br. J. Ophthalmol. 65, 754–759. Fisch, L., Privman, E., Ramot, M., Harel, M., Nir, Y., Kipervasser, S., Andelman, F., Neufeld, M.Y., Kramer, U., Fried, I., Malach, R., 2009. Neural “ignition”: enhanced activation linked to perceptual awareness in human ventral stream visual cortex. Neuron 64, 562–574. Gelbard-Sagiv, H., Mukamel, R., Harel, M., Malach, R., Fried, I., 2008. Internally generated reactivation of single neurons in human hippocampus during free recall. Science 322, 96–101. Goebel, R., 2000. BrainVoyager. Brain Innovation, Masstricht, The Netherlands. Goldwater, 1972. Psychological significance of pupillary movements. Psychol. Bull. 77, 340–355. Golland, Y., Bentin, S., Gelbard, H., Benjamini, Y., Heller, R., Nir, Y., Hasson, U., Malach, R., 2007. Extrinsic and intrinsic systems in the posterior cortex of the human brain revealed during natural sensory stimulation. Cereb. Cortex 17, 766–777. Grandchamp, R., Braboszcz, C., Delorme, A., 2014. Oculometric variations during mind wandering. Front. Psychol. 5, 31. Greicius, M.D., Krasnow, B., Reiss, A.L., Menon, V., 2003. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl. Acad. Sci. U. S. A. 100, 253–258.

D. Yellin et al. / NeuroImage 106 (2015) 414–427 Harmelech, T., Malach, R., 2013. Neurocognitive biases and the patterns of spontaneous correlations in the human cortex. TiCS 10. Hess, E.H., Polt, J.M., 1964. Pupil size in relation to mental activity during simple problemsolving. Science 143, 1190–1192. Hesselmann, G., Kell, C.A., Kleinschmidt, A., 2008. Ongoing activity fluctuations in hMT + bias the perception of coherent visual motion. J. Neurosci. 28, 14481–14485. Hesselmann, G., Flandin, G., Dehaene, S., 2011. Probing the cortical network underlying the psychological refractory period: a combined EEG-fMRI study. NeuroImage 56, 1608–1621. Hurley, L.M., Devilbiss, D.M., Waterhouse, B.D., 2004. A matter of focus: monoaminergic modulation of stimulus coding in mammalian sensory networks. Curr. Opin. Neurobiol. 14, 488–495. Ishai, A., Ungerleider, L.G., Haxby, J.V., 2000. Distributed neural systems for the generation of visual images. Neuron 28, 979–990. Kahneman, D., 1973. Attention and Effort. Prentice-Hall, Englewood Cliffs, N.J. Keller, C.J., Bickel, S., Honey, C.J., Groppe, D.M., Entz, L., Craddock, R.C., Lado, F.A., Kelly, C., Milham, M., Mehta, A.D., 2013. Neurophysiological investigation of spontaneous correlated and anticorrelated fluctuations of the BOLD signal. J. Neurosci. 33, 6333–6342. Kobayashi, M., Imamura, K., Sugai, T., Onoda, N., Yamamoto, M., Komai, S., Watanabe, Y., 2000. Selective suppression of horizontal propagation in rat visual cortex by norepinephrine. Eur. J. Neurosci. 12, 264–272. Kosslyn, S.M., Ganis, G., Thompson, W.L., 2001. Neural foundations of imagery. Nat. Rev. Neurosci. 2, 635–642. Levy, I., Hasson, U., Harel, M., Malach, R., 2004. Functional analysis of the periphery effect in human building related areas. Hum. Brain Mapp. 22, 15–26. Logothetis, N.K., 2010. Neurovascular uncoupling: much ado about nothing. Front. Neuroenerg. 2. Malach, R., Reppas, J.B., Benson, R.R., Kwong, K.K., Jiang, H., Kennedy, W.A., Ledden, P.J., Brady, T.J., Rosen, B.R., Tootell, R.B., 1995. Object-related activity revealed by functional magnetic resonance imaging in human occipital cortex. Proc. Natl. Acad. Sci. U. S. A. 92, 8135–8139. Mukamel, R., Gelbard, H., Arieli, A., Hasson, U., Fried, I., Malach, R., 2005. Coupling between neuronal firing, field potentials, and fMR1 in human auditory cortex. Science 309, 951–954. Murphy, P.R., O'Connell, R.G., O'Sullivan, M., Robertson, I.H., Balsters, J.H., 2014. Pupil diameter covaries with BOLD activity in human locus coeruleus. Hum. Brain Mapp. 35, 4140–4154. Naber, M., Alvarez, G.A., Nakayama, K., 2013. Tracking the allocation of attention using human pupillary oscillations. Front. Psychol. 4. Nir, Y., Mukamel, R., Dinstein, I., Privman, E., Harel, M., Fisch, L., Gelbard-Sagiv, H., Kipervasser, S., Andelman, F., Neufeld, M.Y., Kramer, U., Arieli, A., Fried, I., Malach, R., 2008. Interhemispheric correlations of slow spontaneous neuronal fluctuations revealed in human sensory cortex. Nat. Neurosci. 11, 1100–1108. Quiroga, R.Q., Reddy, L., Kreiman, G., Koch, C., Fried, I., 2005. Invariant visual representation by single neurons in the human brain. Nature 435, 1102–1107. Rajkowski, J., Kubiak, P., Aston-Jones, G., 1994. Locus coeruleus activity in monkey: phasic and tonic changes are associated with altered vigilance. Brain Res. Bull. 35, 607–616.

427

Ramot, M., Wilf, M., Goldberg, H., Weiss, T., Deouell, L.Y., Malach, R., 2011. Coupling between spontaneous (resting state) fMRI fluctuations and human oculo-motor activity. NeuroImage 58, 213–225. Reeves, A., Segal, S.J., 1973. Effects of visual imagery on visual sensitivity and pupil diameter. Percept. Mot. Skills 36, 1091–1098. Rossi, A.F., Paradiso, M.A., 1999. Neural correlates of perceived brightness in the retina, lateral geniculate nucleus, and striate cortex. J. Neurosci. 19, 6145–6156. Salomon, R., Bleich-Cohen, M., Hahamy-Dubossarsky, A., Dinstien, I., Weizman, R., Poyurovsky, M., Kupchik, M., Kotler, M., Hendler, T., Malach, R., 2012. Global functional connectivity deficits in schizophrenia depend on behavioral state. J. Mol. Neurosci. 48, S101–S105. Salomon, R., Levy, D.R., Malach, R., 2013. Deconstructing the default: cortical subdivision of the default mode/intrinsic system during self-related processing. Hum. Brain Mapp. 35, 1491–1502. Sara, S.J., 2009. The locus coeruleus and noradrenergic modulation of cognition. Nat. Rev. Neurosci. 10, 211–223. Schurger, A., Sitt, J.D., Dehaene, S., 2012. An accumulator model for spontaneous neural activity prior to self-initiated movement. Proc. Natl. Acad. Sci. U. S. A. 109, E2904–E2913. Seth, A.K., Chorley, P., Barnett, L.C., 2013. Granger causality analysis of fMRI BOLD signals is invariant to hemodynamic convolution but not downsampling. NeuroImage 65, 540–555. Simpson, H.M., Paivio, A., 1966. Changes in pupil size during an imagery task without motor response involvement. Psychon. Sci. 5, 405–406. Smallwood, J., Brown, K.S., Baird, B., Mrazek, M.D., Franklin, M.S., Schooler, J.W., 2012. Insulation for daydreams: a role for tonic norepinephrine in the facilitation of internally guided thought. PLoS ONE 7, e33706. Soon, C.S., Brass, M., Heinze, H.J., Haynes, J.D., 2008. Unconscious determinants of free decisions in the human brain. Nat. Neurosci. 11, 543–545. Stark, L., Campbell, F.W., Atwood, J., 1958. Pupil unrest: an example of noise in a biological servomechanism. Nature 182, 857–858. Talairach, J., Tournoux, P., 1988. Co-Planar Stereotactic Atlas of the Human Brain. Thieme, Stuttgart, New York. van de Ven, V., Jans, B., Goebel, R., De Weerd, P., 2012. Early human visual cortex encodes surface brightness induced by dynamic context. J. Cogn. Neurosci. 24, 367–377. Wang, X.J., 2008. Decision making in recurrent neuronal circuits. Neuron 60, 215–234. Wang, L., Hermens, D.F., Hickie, I.B., Lagopoulos, J., 2012. A systematic review of restingstate functional-MRI studies in major depression. J. Affect. Disord. 142, 6–12. Weissenbacher, A., Kasess, C., Gerstl, F., Lanzenberger, R., Moser, E., Windischberger, C., 2009. Correlations and anticorrelations in resting-state functional connectivity MRI: a quantitative comparison of preprocessing strategies. NeuroImage 47, 1408–1416. Wierda, S.M., van Rijn, H., Taatgen, N.A., Martens, S., 2012. Pupil dilation deconvolution reveals the dynamics of attention at high temporal resolution. Proc. Natl. Acad. Sci. U. S. A. 109, 8456–8460. Zylberberg, A., Oliva, M., Sigman, M., 2012. Pupil dilation: a fingerprint of temporal selection during the “attentional blink”. Front. Psychol. 3, 316.