Spontaneous eye-movements during eyes-open rest reduce resting-state-network modularity by increasing visual-sensorimotor connectivity

During wakeful rest, individuals make small eye movements when asked to fixate. We examined how these endogenously-driven oculomotor patterns impact topography and topology of functional brain networks. We used a dataset consisting of eyes-open resting-state (RS) fMRI data with simultaneous eye-tracking (Nilsonne et al., 2016). The eye-tracking data indicated minor movements during rest, on the order of 1.0 degree on average when analyzed over 2sec epochs, which correlated modestly with RS BOLD data. However, the eye-tracking data correlated well with echo-planar imaging (EPI) time series sampled from the area of the Eye-Orbit (EO-EPI), which is a signal previously used to identify eye movements during exogenous saccades and movie viewing. We found that EO-EPI data correlated with activity in an extensive motor and sensory-motor network, but also some components of the dorsal attention network including the frontal and supplementary eye fields. Partialling out variance related to EO-EPI from RS data reduced connectivity, primarily between sensory-motor and visual areas. For three different network sparsity levels, the resulting RS connectivity networks showed higher modularity, lower mean connectivity strength, and lower mean clustering coefficient. Our results highlight new aspects of endogenous eye movement control during wakeful rest. They show that oculomotor-related contributions form an important component of RS network topology, and that those should be considered in interpreting differences in network structure between populations, or as a function of different experimental conditions.


14
The study of human brain activity during resting state (RS) is of considerable interest in both basic and clinical brain research.  rest. For this reason, researchers often remove effects of motion and physiology from RS data prior to analysis. 30 Here we examined how RS connectivity is related to a different factor, which is eye movement during rest. For purposes 31 of understanding endogenous computations, spontaneous eye-movement at rest straddles the space between an interesting 32 neurobiological phenomenon reflecting the output of endogenous activity and a nuisance factor reflecting motor activity. On 33 Specific aims (1) f 2 = 1/(pupilwidth 2 ) + 1/(pupilheight 2 ) (2) After transforming the using these functions, we defined artifacts as sharp peaks in the resulting time series. In addition, we 125 defined blinks as artifacts whose duration was 100-400 msec. After removing time points containing artifacts, we considered 126 the time series of gaze locations. To limit the influence of the noise due to data acquisition failures we only included data 127 from participants with gaze variance below an arbitrary value of 0.1 radians. Consequently, we analyzed data from 32 (of 128 77) participants. For these, the proportion of artifacts was on average 18 ± 2%; blinks occurred with an average period of First, we applied brain extraction and tissue segmentation (Gray Matter, White Matter, CSF) to the structural T1 images using the antsBrainExtraction function of ANTs software (Avants, Tustison, & Song, 2011). We used ANTs for all registration 134 routines in our pipeline. We registered each participant's structural image to standard space using non-linear registration (ICBM 135 2009 non-linear assymetric template; Fonov, Evans, McKinstry, Almli, & Collins, 2009), and saved the inverse of the warps. 136 We also registered the structural and functional images using affine transformation. We used the combination of these two 137 transformations to align data from each participant's original space to common space, or vice versa, in a single step. 138 To delineate each participants "eye orbit" area we first marked this area on the common-space template. We then transformed 139 this mask to each participant's original space, and made any additional modifications (if needed) therein. Specifically, we 140 delineated anatomical masks of the "eye orbit" area in common space using MRICRON (Rorden, Karnath, & Bonilha, 2007), 141 for which we used an MNI template provided with FSL (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012). Both 142 eye orbits were included in the mask. The masks' location was transformed to each participant's individual space using 143 the combination of T1→subject space matrices and inverse of the T1→MNI matrices mentioned above. We also created 144 cerebral-spinal fluid (CSF) and white matter masks in MNI space and transformed them to individual space, where they were 145 eroded by one voxel from each direction to be more conservative. We then extracted the mean time series from these white 146 matter and CSF masks. These were used as nuisance regressors in an initial regression (details below). 147 We used AFNI (Cox, 1996) for analyzing the functional RS images. We implemented the following steps: slice time 148 correction, motion correction (base image: first volume of the run), and band-pass filtering (0.01 − 0.1Hz). To remove other 149 nuisance sources of variance from the functional times series we implemented preliminary data-cleaning using regression with 150 the following regressors: i) motion parameters estimated during motion correction, ii) mean white matter and CSF time series,

151
iii) and frame-wise displacement values (included in the model as a regressor). We considered the residuals of this regression as 152 a "cleaned" time series that was the starting point for further analyses.

153
To improve signal to noise of the subsequent regression models which were of primary interest, we then spatially-smoothed 154 the cleaned time series with a 6mm FWHM kernel. From this time series we also derived an Eye-Orbit EPI regressor, which was 155 defined as the mean time series from both eye-orbit masks, after spatial smoothing, which we refer to as EYE raw . We convolved 156 the EYE raw with an HRF basis function (Using AFNI's waver command), producing a EYE conv time series. In separate analyses 157 we used either EYE raw or EYE conv as "seed" regressors, to identify EO-EPI-correlated brain areas.

158
Determining correlation between eye-tracking measures and EO-EPI time series 159 We were interested in the relationship between several measures of eye movement and the EPI time series sampled from the 160 eye-orbit regions (EO-EPI series). We derived 12 time series from the eye-tracking data were: the measured gaze location,

161
GazeX these correlation values at the group level using a T-test.

Correlates of Eye-tracking metrics
We examined whole-brain RS correlations with several eye tracking measures: GazeX, GazeX 2 , vel GazeX, vel GazeX 2 , Pupil size, and blink function. The BOLD data modeled were "cleaned" time series from which only typical artifact sources 189 were removed. We implemented two modeling approaches: In the first, we resampled each eye-tracking measure of interest to 190 the sampling resolution of the MR acquisition (Hz = 0.4) and convolved the result with canonical HRF via AFNI's Waver 191 function to construct a regressor. In the second, we used a Finite Impulse Response (FIR) function modeling approach where 192 the BOLD impulse response was estimated using six tent functions (using AFNI's tent basis function). This approach does not 193 assume a fixed shape. From these estimates, we averaged the first three beta coefficients (corresponding to 0 − 7.5sec post 194 eye-tracker dynamics) and propagated the value to a group-level analysis. Family wise error correction was implemented using 195 FSL's TFCE implementation.

197
Beta values associated with EYE conv or EYE raw were transformed to MNI space. To identify clusters where these beta values 198 were significantly positive or significantly negative we computed voxel-wise statistics (Wilcoxon signed-rank test) and then 199 implemented cluster-level control for family-wise-error using permtuations as described below. We used a non-parametric test were either positive or negative). We saved the largest cluster size from each permutation, and the resulting set of 10,000 values 209 of largest-cluster sizes defined the sampling distribution. The 95% percentile rank entry of the sampling distribution served as 210 the critical value. This value was used to define statistically-significant clusters in the experimental data.

211
To evaluate whether significant EO-EPI correlates were found in areas dominated by artifacts, we calculated temporal signal 212 to noise ratio (tSNR) for each participant. To create tSNR map for each participant, we used the raw functional images (before 213 applying any signal processing steps), but after removal of the initial 10 stabiliziation images. We divided the absolute mean 214 value of each voxel by its standard deviation. We then applied the statistically significant clusters found for EYE raw and EYE conv 215 series as masks to determine mean of the tSNR in each statistically significant spatial cluster. The motivation for this analysis 216 was a prior report (W. Chen & Zhu, 1997) showing that Nyquist ghosting artifacts can propagate eye signals into midbrain 217 areas (in the case of axial acquisition). Two MR physicists examined the QA reports produced by the scanner and did not find 218 evidence for ghosting. However, we still wanted to evaluate if any EO-EPI whole-brain correlates were found in regions with 219 low tSNR. those two ROIs extracted from the cleaned and smoothed data. We constructed a regression model to predict that FEF and SEF 227 regional activity from the EO-EPI series per participant. Coefficients were analyzed using a Wilcoxon rank sum test. where they were further limited to gray matter by multiplying all ROIs with the participant-specific gray matter mask (to limit 232 the influence of data from non-gray matter areas). We extracted the mean time series from each ROI, for the two types of 233 spatially-smoothed resting-state data we derived (one typical, and the other with EO-EPI regressed). We examined the network 234 features after thresholding the connectivity matrices at three sparsity levels -10%, 20%, 30%. While the thresholding removed 235 weak connections, we did not binarize the remaining (above-threshold) connections but analyzed the complete set of data.

236
From each participant's resting state network we derived the following metrics: node degree, strength, cluster coefficient, 237 transitivity, assortativity, efficiency, number of communities, betweenness centrality and modularity. We calculated these using 238 the Brain Connectivity Toolbox (Rubinov & Sporns, 2010) (See Appendix for description of the metrics as described in the 239 Brain Connectivity Toolbox). We calculated these parameters for the original and "clean" networks as defined above. We then 240 tested which of these parameters differed as a result of the EO-EPI-removal procedure using paired-sample T-tests. We defined 241 a robust result as one that was statistically significant across all three levels of sparsity thresholding (full results are reported for 242 completeness). 243 We also probed changes in global topology by quantifying the impact of EO-EPI removal on the shape of the entire degree thresholds. These criteria required that the value of a node be higher than 1 SD above the mean value for each of these 251 empirical distributions: node strength, node degree and node between-ness centrality. Nodes matching all three criteria were 252 considered hubs. The chance probability of a node being a hub (assuming a normal distribution) is ∼ 0.34 3 = .04. To evaluate 253 whether removal of EO-EPI variance impacted whether a region satisfied hub criteria, for each region we counted the number 254 of participants for which the region was classified as a hub, with our without EO-EPI removal. On a binomial, a difference 255 would need to consist of at least 7 or more participants (binomial test parameters: N = 86; K = 7; p = .04). 256 We also identified any specific pair-wise differences in regional connectivity for the raw and cleaned matrices. After 257 applying Fisher's Z transformation, pair-wise correlation values were subjected to paired-sample t-tests. We used false discovery 258 rate (FDR) to correct for multiple comparisons.

259
Dual Regression 260 We used dual regression to determine how removal of activity associated with the EO-EPI regressor impacted connectivity in 261 well-defined resting-state networks. The procedure was implemented in AFNI and followed workflows described previously removed, and one from the functional data from which this variance was removed using the EYE conv regressor.

268
To determine whole brain connectivity of the seed regions we inserted all 14 time series into a single multiple regression.

269
In effect, we conducted two separate regression models: Model #1 was a "typical" model where the mask-derived seed time 270 series produced from the original (typically-processed) functional data served as regressors to predict whole-brain resting state 271 data. This process replicates the procedure typically used in resting-state dual regression. Model #2 was an "EO-EPI-removed" 272 model where masked-derived seed time series produced from EO-EPI-removed functional data were used as regressors to 273 predict the EO-EPI-removed functional data. This is a dual regression that uses EO-EPI-cleaned resting state data.

274
The produced beta weights were analyzed using group level repeated-measures test to identify seed-time-series whose 275 connectivity differed between the two data sets; i.e., whose connectivity was impacted by the EO-EPI-removal procedure. We 276 used FSL's randomise function (Jenkinson et al., 2012). A within group T-test with 10000 permutations and threshold-free 277 cluster enhancement was applied. Because our interest was in evaluating the impact EO-EPI-regressor we adopted a liberal 278 approach of not correcting for multiple comparisons across the 14 networks tested in the dual regression procedure. We also 279 note that the 14 time series used for dual regression were relatively weakly correlated in this data set: to determine collinearity, 280 on the single participant level we computed the 14 × 14 cross-correlation matrix and then averaged these across participants.

281
The highest mean correlation was 0.55, which licensed separate analyses for each network regressor.

283
Eye tracking data: Quality and correlation with whole-brain BOLD 284 Based on our artifact rejection criteria, usable eye-tracking data were available for 32 of 77 participants for which eye tracking Each analysis is corrected for multiple comparisons using FSL's implementation of TFCE Family-wise-error control.
in Figure 1B regressor. Regressions based on canonical HRF-convolved regressors produced results that were not statistically significant.

295
Eye tracking data: Correlation with Eye Orbit EPI data 296 We evaluated the correlation between each of the 12 types of eye tracking time series (see Methods) and the EO-EPI data. We 297 controlled for the 12 tests using Bonferroni correction, because some of the tests were highly correlated (see Figure Appx.2). 298 We found that three eye-tracking regressors significantly correlated with the EO-EPI envelope (Bonferroni corrected for 12 299 tests): the gaze power vel GazeX 2 + vel GazeY 2 , square of pupil size PupilSize 2 , and the gaze velocity in the vertical (Y) 300 direction. The pupil size was evaluated as deviation from the subject's mean value, so its squared value indicated absolute 301 deviations from mean value (we used squared deviation rather than absolute value as the derivative of the exponent is better 302 behaved than that of the absolute function). Figure 2A shows sample time series reflecting raw EO-EPI, its envelope and 303 eye-tracking regressors, and Figure 2B shows the estimated Kernels for gaze power and square of pupil size.

304
Pupil-size squared explained 7 ± 2% of the variance of the EO-EPI envelope and presented a significant positive correlation:

Connectivity of EO-EPI regressors
We identified an extensive system that correlated with the EO-EPI regressor. For the convolved version of the EO-EPI regressor 314 (EYE conv ) we found correlations in pre-and post-central gyri bilaterally, parts of the superior temporal gyrus and visual cortex 315 ( Figure 3A). We also identified strong correlations (of opposite sign) in the thalamus ( Figure 4A). We also found whole brain  Supplementary Table 5). For EYE conv the clusters found in the left and right cerebellum 326 were associated with low tSNR (and relatively systematically across participants, see Supplementary Table 6), as was a cluster 327 in the mid occipital gyrus bilaterally (potentially as it includes time series from the field of view between the two hemisphere).  The reason for these differences across participants is unclear. However, a byproduct is that when the EO-EPI regressor

344
The large standard deviation of the EO-EPI regressor was related to peaks in that signal. As indicated in the Methods 345 section, applying a 'despiking' procedure reduced the sensitivity of the whole brain correlation analysis: its most extreme   Fitting the degree distributions using an exponentially truncated power law showed that the EO-EPI removed networks 366 differed in the degree distribution (see Figure 6). As shown in the Figure  No statistically significant differences were found for 30% sparsity networks. Figure Appx.6 presents mean degree-distributions 371 for Raw and Clean networks in the different sparsity levels. 372 We determined which areas tended to show changes in connectivity as function of EO-EPI removal. In general, this analysis 373 is not independent of the whole-brain correlation with the EO-EPI time series used as a regressor, but it is more sensitive  which showed stronger connectivity with multiple other brain areas.

380
The dual regression analysis did not identify any pre-defined RS network for which connectivity changed significantly.

381
A hub-focused analysis that examined whether there were regions more frequently identified as hubs in the raw or EO-EPI-382 removed series also produced a null result: the most extreme example was a region defined as hub for 20 participants in one case 383 and 25 in another (a non-significant difference on a binomial). While the location of these hubs was not a central point of the 384 current study, broadly speaking, for the 10% sparsity threshold (raw) matrices, hubs were localized motor and sensory-motor 385 areas (9 regions) Dorsal attention (6 regions), DMN (4 regions), temporal-parietal areas (4 regions) and ventral attention (2 386 areas). Only one visual extrastriate area was identified as a hub.

388
Neuroimaging is continuously expanding our understanding of the principles that determine structured patterns of RS connec-389 tivity. Our findings demonstrate that endogenous eye movements during RS contribute significantly to structured patterns of RS 390 connectivity. Our main finding is that eye movements, measured via EPI time series recorded from the eye orbits, identified a 391 sensory-motor system that appeared linked to oculomotor activity. Removal of activity accounted for by eye movements had 392 systematic impact on whole-brain connectivity. We first address issues related to oculomotor measurement during the resting 393 state that emerged in the study and then discuss the implications of the results for basic and applied research.

394
Probing resting-state networks with Eye tracking and eye-orbit EPI data: technical considerations 395 As reviewed in the Introduction, few studies have studied brain activity patterns that are correlated with oculomotor activity data, which presented itself in higher power across all frequencies for rejected data as compared to analyzed data. 402 We found correlations between the eye-tracking metrics and EPI data recorded from the eye orbit area (EO-EPI), Bonferroni 403 corrected for 12 correlation tests. These were found for Gaze power, pupil size (squared), and gaze velocity in the Y direction.

404
These data are consistent several prior reports. Beauchamp (2003) showed that peaks in the EO-EPI time series occur when an  coupling in which changes in pupil size preceded EO-EPI fluctuations (the latter delayed by 2 sec), and of a strong peak 416 frequency of 0.04Hz for EO-EPI are both consistent with the possibility that EO-EPI also reflects metabolic activity. We also 417 found little independent evidence to suggest a strong contribution of motion artifacts to EO-EPI: beyond one participant for 418 which 5 of 6 motion parameters correlated with EO-EPI, we only found 2 additional correlations with motion elements, for two 419 additional participants.

420
Note that task compliance during this RS study was good. First, participants were continuously monitored and experimenters 421 verified participants did not drift off to sleep during the scan. Second, the eye-tracking data indicated compliance with the 422 task instructions in that the eye movements that were made during fixation were minor in magnitude (see Figure 1A). When 423 evaluating average movements between successive 2sec epochs we found that in 75% of the cases, the magnitude was below 424 2degree, which corresponds to a small displacement. For this reason, we consider these data to be representative of typical 425 compliant behavior during wakeful rest.

427
When used as a whole-brain regressor, the EO-EPI time series correlated with an extensive bilateral sensory-motor system. In    Figure Appx.1. Power spectra of eye-tracking data for data rejected or maintained by the quality-control procedure. Note the dual Y-axis scales. Rejected data presented power around one order of magnitude higher than maintained data. This difference was significant even for the highest measurable frequency, where power in rejected data was 2.5 times higher than that maintained.

521
• Degree: Node degree is the number of links connected to the node. In directed networks, the in-degree is the number of 522 inward links and the out-degree is the number of outward links. Connection weights are ignored in calculations.

523
• Strength: Node strength is the sum of weights of links connected to the node. In directed networks, the in-strength is the 524 sum of inward link weights and the out-strength is the sum of outward link weights.

525
• Clustering coefficient: The clustering coefficient is the fraction of triangles around a node and is equivalent to the 526 fraction of node's neighbors that are neighbors of each other.

527
• Transitivity: The transitivity is the ratio of triangles to triplets in the network and is an alternative to the clustering 528 coefficient.

529
• Assortativity: The assortativity coefficient is a correlation coefficient between the degrees of all nodes on two opposite 530 ends of a link. A positive assortativity coefficient indicates that nodes tend to link to other nodes with the same or similar 531 degree.

532
• Global Efficiency : The global efficiency is the average inverse shortest path length in the network, and is inversely 533 related to the characteristic path length.

534
• Number of communities: The optimal community structure is a subdivision of the network into nonoverlapping groups 535 of nodes in a way that maximizes the number of within-group edges, and minimizes the number of between-group edges.

536
To calculate maximum number of community for each network, we recorded the largest number of community.

537
• Modularity: The modularity is a statistic that quantifies the degree to which the network may be subdivided into such 538 clearly delineated groups.

539
• Betweenness Centrality: Node betweenness centrality is the fraction of all shortest paths in the network that contains a 1 This mis-measurement is known as the Pupil Foreshortening Error (PFE). For example, Hayes and Petrov (2016) showed that deviations from center of 547 camera-view produce systematic PFEs that can reach ∼ 12% at typical viewing distances. Significant PFEs were produced even with movements as small as 548 ∼ 4 • of center. 549 2 All data available online at https://www.openneuro.org/datasets/ds000201/ 550 3 Pupil size was defined as (pupil width + pupil height)/2; we note that with our instrumentation (as well as many other eye trackers), pupil size is 551 confounded with gaze position (Hayes & Petrov, 2016), resulting in significant correlations between pupil size and the gaze both in x and y directions (p < .01 552 in 30/32 subjects).