Anticorrelated inter-network electrophysiological activity varies dynamically with attentional performance and behavioral states

The default mode network (DMN) is thought to exhibit infraslow anticorrelated activity with dorsal attention (DAN) and salience (SN) networks across various behavioral states. To investigate the dynamics of activity across these networks on a finer timescale, we used human intracranial electroencephalography with simultaneous recordings within core nodes of the three networks. During attentional task performance, the three sites showed dissociable profiles of high-frequency broadband activity. Anticorrelated infraslow fluctuations of this activity were found across networks during task performance but also intermittently emerged during rest and sleep in concert with the expression of task-like network-level topographic patterns. Critically, on a finer timescale, DAN and SN activations preceded DMN deactivations by hundreds of milliseconds. Moreover, greater lagged, but not zero-lag, anticorrelation between DAN and DMN activity was associated with better attentional performance. These findings have implications for interpreting antagonistic network relationships and confirm the behavioral importance of time-lagged inter-network interactions.


Introduction
The brain structures important for attention have long been described as a set of discrete components of networks, each of which serves a unique role in cognitive interactions with the sensory environment (Posner and Petersen, 1990). The dorsal attention network (DAN) is implicated in goal-directed (top-down) attention (Corbetta and Shulman, 2002), and the salience network (SN) is implicated in stimulus-driven (bottom-up) attention and cognitive control (Seeley et al., 2007, Uddin, 2015. Given functional neuroimaging evidence, these described networks exhibit opposing activity with yet another set of regions that constitute the default mode network (DMN) (Raichle et al., 2001, Fox et al., 2005, which unlike DAN and SN structures, tends to show deactivation during conditions involving externally-oriented attention (Buckner et al., 2008).
It has been hypothesized that anticorrelated activity between the DMN and DAN/SN reflects a competition for control over shared computational resources and could be a marker of one's level of engagement with externally-oriented tasks (Sonuga-Barke and Castellanos, 2007).
In addition to task-dependent activity, fMRI studies during wakeful rest have shown spontaneous anticorrelated DMN-DAN/SN activity in infraslow (<0.1 Hz) fluctuations of blood-oxygen-leveldependent (BOLD) signals (Fox et al., 2005, Fransson, 2005) -a finding that has remained contentious (Murphy and Fox, 2017). Persistence of such anticorrelated activity in task-free states would potentially suggest that functionally competing systems, characterized by continual switching between internally-and externally-biased modes of attention, are an intrinsic property of the brain (Buckner et al., 2013, Honey et al., 2018. To date, however, evidence for anticorrelated activity has relied almost exclusively on fMRI, which offers limited temporal resolution and necessitates data preprocessing that may bias estimates of negative correlations Glover, 2009, Murphy et al., 2009). Intracranial electroencephalography (iEEG) in human subjects offers anatomical precision, high temporal resolution, and sensitivity to activity in the high-frequency broadband (HFB, also known as high gamma) range (~70-170 Hz)-a well-established correlate of the BOLD signal and neuronal population spiking (Parvizi and Kastner, 2018). A handful of iEEG studies involving recordings from putative DMN and DAN nodes have shown task-evoked HFB responses that resemble the antagonistic inter-network patterns observed in fMRI (Ossandon et al., 2011, Ramot et al., 2012, Raccah et al., 2018. In addition, resting state iEEG has revealed correlates of the DMN and DAN (Foster et al., 2015, Hacker et al., 2017, Kucyi et al., 2018a and that a subset of region pairs with resting BOLD anticorrelations exhibit weaker, but significant anticorrelations of slow (0.1-1 Hz) HFB activity (Keller et al., 2013).
The extant iEEG evidence provides promising initial electrophysiological validation of antagonistic inter-network interactions, but beyond confirming fMRI findings, critical open questions remain: Is infra-slow anticorrelated activity between task-responsive DMN, DAN and SN neuronal populations found in task-free states and in specific frequency components of electrophysiological signals? Do nodes of the DMN, DAN and SN exhibit distinguishable, systematic temporal profiles of task-evoked electrophysiological activity? Are time-resolved antagonistic interactions relevant to intra-individual variations in attentional task performance?
Here we report an iEEG investigation of anticorrelated brain networks across multiple sessions of attentional task performance, wakeful rest and sleep. We localized iEEG recording sites in key nodes of the DMN, DAN, and SN through a rigorous survey of anatomical boundaries, iEEG response profiles, and within-individual resting-state fMRI connectivity. We then examined how inter-network iEEG interactions vary as a function of attentional (task performance) and behavioral (active-task versus rest versus sleep) states. We hypothesized that temporal features of electrophysiological anticorrelated activity would be associated with fluctuations in sustained attention and that the overall magnitude of anticorrelated activity would co-vary within and across behavioral states despite a constrained wider network-level spatial organization.

Unique iEEG cohort
We obtained iEEG recordings from a total of 896 recording sites in seven participants (S1-S7, five with depth electrode and two with subdural recordings) undergoing treatment for focal epilepsy.
Data reported here are from regions void of pathological activity and outside each participant's seizure zone. An average of 78±14% of channels per participant were retained for analysis, a percentage that is within the typical range (Parvizi and Kastner, 2018). Each subject had simultaneous electrode coverage across key nodes of DMN, DAN and SN, as defined based on anatomical boundaries and, when available (n=5), confirmed with individual-level resting-state fMRI (see Methods). Specifically, we focused on recordings from three regions: 1) posteromedial cortex (PMC) within the DMN (Raichle et al., 2001); 2) dorsal posterior parietal cortex (dPPC), including superior parietal lobule and intraparietal sulcus, within the DAN (Corbetta and Shulman, 2002); and 3) dorsal anterior insular cortex (dAIC) within the SN (Seeley et al., 2007). Our cohort was unique and carefully selected from a larger group of patients based on the presence of recording sites within the three regions of interest. Implantation of intracranial electrodes in all subjects were solely based on clinical needs, and the location and number of electrodes implanted in each case were decided by a collective agreement of medical staff at Stanford Medical Center (6 subjects) and Beijing Tiantan Hospital, Capital Medical University (1 subject).
Subjects performed between four to eight sessions (total duration range: 24-48 minutes per subject) of the Gradual-onset Continuous Performance Task (GradCPT), a test of sustained attention that has reliably been associated with anticorrelated DMN versus DAN/SN activity in fMRI studies (Esterman et al., 2013, Kucyi et al., 2016, Fortenbaugh et al., 2018. Gradually changing images of scenes were presented every 800 ms, and subjects were instructed to respond with a button press to city (frequent) but not to mountain (infrequent) scenes ( Figure   1A). The task requires sustained attention and withholding of a habitual response to infrequent target events. We defined behavioral performance within each session of the experiment by the measure of sensitivity (d′), which is based on the accuracy of task performance (accounting for both hits and false alarms) (Fortenbaugh et al., 2015, Rosenberg et al., 2016. We found that performance varied from session to session within subjects (Table S1). In three subjects, additional recordings were obtained during wakeful rest and sleep sessions of similar durations to those obtained during task performance.

Functional localization of electrophysiological antagonistic networks
We functionally localized electrode contacts of interest within the PMC, dPPC and dAIC anatomical boundaries ( Figure 1B, see Methods for details). These sites were chosen as core nodes of three networks of interest, namely the DMN, DAN, and SN respectively. We analyzed electrophysiological activity during withheld presses to infrequent targets (correct omissions) relative to presses for frequent non-targets (correct commissions) ( Figure 1C).
In our analysis, we adopted a recently described data-driven multivariate approach to decompose the relative contributions of distinct bands of the iEEG signal to task-evoked responses (Schrouff et al., 2016). Specifically, for correct omission (mountain) and correct commission (city) trials, we generated kernels of similarity across trials for each channel that was anatomically within at least one of the ROIs (dPPC, PMC, dAIC) and the power amplitudes of seven frequency bands [δ (1-3 Hz), θ (4-7 Hz), α (8-12 Hz), β1 (13-29 Hz), β2 (30-39 Hz), γ (40-70 Hz) HFB (70-170 Hz)]. Using these features (i.e., 7 frequencies x number of ROI channels in a given subject), multikernel learning analyses revealed significant classification accuracy within each subject (M±SD=80.7±7.0%, p<0.01 in all cases) ( Figure 1E). Across frequencies, HFB features had the highest mean contributions to classification accuracy across subjects, although lower frequencies also contributed ( Figure 1F). Multikernel learning models based on single ROIs (and the seven frequencies) revealed that classification accuracy was strongest for dPPC (M=79.3%) and dAIC (M=78.5%) and weakest for PMC (M=66.1%) activity ( Figure 1G). Based on these findings, we focus our central further iEEG analyses on HFB activity but also consider other frequency ranges.
As noted, HFB (~70-170 Hz) is a well-established correlate of the BOLD signal (and neuronal population spiking) (Parvizi and Kastner, 2018). Based on previous fMRI findings (Fortenbaugh et al., 2018), we expected decreased HFB activity in the PMC following targets relative to nontargets, while in dPPC and dAIC, we expected greater HFB activity following targets relative to non-targets. In each subject (n=7), we successfully identified peak-responsive channels within the PMC that showed temporal clusters of significantly decreased HFB power following target compared to non-target onsets (Monte Carlo p<0.002 in all instances, corrected for number of channels within ROI using cluster-based permutation testing). In each participant with dPPC (n=6) and dAIC (n=4) coverage, we also identified temporal clusters of significantly increased HFB power within the dPPC (Monte Carlo p<0.001 in all instances) and dAIC (Monte Carlo p<0.001 in all instances) for target compared to non-target trials (see Figure S1A for all subjects). In the peak-responsive channels, spectrograms suggested that deactivations (PMC) and activations (dPPC, dAIC) were consistent within the HFB range, although responses were also evident in lower frequency ranges (See Figure S1B for all subjects). In addition to demonstrating task-evoked iEEG response profiles, we confirmed the network membership of peak-responsive channels within each ROI using resting-state fMRI. To do so, we extracted the BOLD time series from the channels' locations and performed seed-based functional connectivity analysis (Biswal et al., 1995) to identify remote brain regions with correlated infraslow BOLD activity. In the five subjects who underwent resting-state fMRI, wellknown features of the DMN, DAN and SN were found at the individual level (single-subject example in Figure 1D; see Figure S1C for all subjects). Specifically, for the PMC peak-deactive iEEG sites, selectively correlated BOLD activity was found with established nodes of the DMN including medial prefrontal cortex and angular gyrus. For the dPPC peak-active iEEG sites, correlated BOLD activity was found with nodes of the DAN including frontal eye fields and area MT+. For dAIC peak-active iEEG sites, correlated BOLD activity was found with the mid-cingulate cortex and anterior temporoparietal junction within SN. Also apparent in these BOLD network maps were DMN-DAN and DMN-SN anticorrelations, which were obtained using global signal regression but which also were apparent when using an alternative preprocessing pipeline ( Figure S2).
Overlays of DMN, DAN and SN templates, derived from a population-level atlas (Yeo et al., 2011), showed concordance between BOLD networks in each individual subject and those found in neurotypical adults (black outlines on cortical surfaces in Figures 1D and S1C). In the two subjects that did not undergo fMRI, registration of the Yeo atlas to individual cortical surfaces suggested that peak-responsive PMC, dPPC, and dAIC channels were respectively within the boundaries of the fMRI-defined DMN, DAN and SN. The correspondence between iEEG response profile and fMRI network identity could be illustrated in cases where subdural electrodes densely covered areas near network boundaries. For example, in S5, dPPC and PMC responsive sites were found within the Yeo atlas' DAN and DMN, respectively, but HFB responses at adjacent sites a few millimeters away were absent or distinct ( Figure 1H). These network-wide spatial profiles, together with the observed task-evoked iEEG response profiles, suggest that the general functions and connectivity of the PMC, dPPC and dAIC were preserved within participants.

Figure 1. Task paradigm and functional localization of electrode sites in three regions of interest.
A) The Gradual-Onset Continuous Performance Task. City and mountain scene images faded continuously from one image to the next every 800 msec. Subjects were instructed to press a button when they saw cities but to withhold responses to (rare) mountains. Trial onset (orange arrow) was the time at which stimulus fade-in was initiated. B) Anatomical locations of electrode contacts implanted in three regions of interest that showed peak-responsive high-frequency broadband (HFB; Figure S1 and Figure S2.

Dissociable DAN and SN response profiles
As noted, both dPPC and dAIC exhibited similarly increased electrophysiological responses in the HFB range during correct omission trials ( Figure S3A). However, during commission errors that signify lapses of attention (i.e., when participants incorrectly pressed the key button to infrequent mountain targets), we observed a clear divergence between DAN and SN sites. As expected from fMRI studies (Ham et al., 2013, Neta et al., 2015, Fortenbaugh et al., 2018, during commission errors compared with correct omission trials, the dAIC showed temporal clusters of significantly greater HFB activation in all 4 subjects with dAIC coverage (Monte Carlo p<0.05 in all cases) whereas the dPPC showed no significant differences in any of the 6 subjects with dPPC coverage (Figure 2A and 2C). The error-related HFB increases in the dAIC appeared later and more sustained compared with those found for correct omission trials ( Figure 2C and 2D). The PMC showed largely similar HFB deactivation for correct omission and commission error trials, although in 2 out of 7 PMC subjects, temporal clusters of greater deactivation for correct trials were found (Monte Carlo p<0.05 in both cases) ( Figure 2B). Taken together, these findings further confirm that dPPC and dAIC recording sites were within functionally dissociable networks. . D) Location of peak-responsive dAIC electrode and spectrograms for correct and error trials in an example subject (S1). For error compared to correct trials, a stronger and more sustained increase in power amplitude can be seen in the HFB range.

Diminished anticorrelated iEEG activity during rest and sleep states
Past fMRI studies have suggested the possibility of anti-correlated inter-network activity during non-task states (Fox et al., 2005, Fransson, 2005 -but these findings have remained disputed (Murphy and Fox, 2017). Since task-independent anticorrelated inter-network activity in fMRI has relied on infraslow (<0.1 Hz) fluctuations of BOLD signals, we used comparable parameters to explore the presence of task-independent intrinsic anticorrelations between PMC and dPPC and between PMC and dAIC. For this, we selected the same functionally localized (peak-responsive) PMC, dPPC and dAIC sites shown in the previous analyses and assessed whether infraslow (<0.1 Hz) HFB envelope anticorrelation was found between them.
In three patients (S1, S2, S6) with rest and sleep ( Figure S4)  independent 100-second windows and compared infraslow HFB functional connectivity across states. We found that within each subject, there was a significant interaction between behavioral state and infraslow HFB functional connectivity for dPPC-PMC (F>10, p<0.001 in all cases) and dAIC-PMC (F>10, p<0.001 in all cases) channel pairs. Though some time windows with anticorrelated activity were found during rest and sleep, consistent anticorrelation across windows was reliably detected only in the task state for both dPPC-PMC ( Figure 3A) and dAIC-PMC ( Figure 3C) channel pairs. In all subjects, infraslow HFB functional connectivity in task was significantly lower than that in rest for dPPC-PMC (p<0.005 in all cases) and dAIC-PMC (p<0.005 in all cases) channel pairs. Infraslow HFB functional connectivity was also significantly lower in task than in sleep for dPPC-PMC (p<0.001 in all cases) and dAIC-PMC (p<0.001 in all cases) channel pairs. Differences between rest and sleep were found less consistently for both dPPC-PMC (S1: p=0.03; S2: p=0.14; S6: p=0.81) and dAIC-PMC (S1: p=0.89; S2: p=0.046) channel pairs.
To assess whether infraslow anticorrelated activity varied as a function of carrier frequency, we compared HFB with lower frequency ranges within subjects. We found that task anticorrelation was strongest and most consistent in the HFB range for both dPPC-PMC ( Figure 3B) and dAIC-PMC ( Figure 3D) channel pairs. Lower frequency ranges instead showed positive, rather than negative, correlations between regions across task, rest and sleep states.

Task-like topographic network patterns during spontaneous iEEG anticorrelations
However, network-level topographic patterns also show spontaneous changes across various time scales, suggesting that the brain continually shifts across various 'states,' and some of the states could show more task-like activity patterns than others (Deco et al., 2013). As we found that the presence and magnitude of anticorrelated activity was variable across time windows in rest and sleep, we hypothesized that network-level topographic patterns would be most similar to task-like patterns when spontaneous anticorrelated activity was increased. This could suggest that spontaneous anticorrelations signify the emergence of vigilant 'task-like' states.
We first confirmed that a relatively stable spatial topography of activity interactions was present across all states. To do so, we performed a wider network-level analysis of correlated infraslow HFB activity from functionally localized dPPC, PMC and dAIC channels (seed regions) to all other implanted electrode contacts (target regions) within a given subject's brain (excluding channels that were immediately neighboring the seed, were deemed pathological, were in white matter, or percentile; see highlighted data points in Figure 4).  Figure S5). Thus, our results are consistent with the notion that intrinsic functional connectivity remains relatively stable across states despite the presence of cross-state changes in the magnitude of inter-network anticorrelation.
We next tested whether spontaneous temporal fluctuations in the degree of inter-network anticorrelation was associated with variation in the 'task-like' quality of topographic network patterns. In each 100-second window in rest and sleep, we determined the magnitude of dPPC-PMC or dAIC-PMC infraslow HFB anticorrelation, and we extracted the all-to-all connectivity matrix for implanted electrodes. We then compared the similarity of each window's matrix with a 'task template' matrix that was constructed based on the average all-to-all connectivity pattern during continuous task performance ( Figure 5A). This analysis revealed that temporal windows at rest with greater rest-task topographic similarity were associated with greater resting state anticorrelation between dPPC and PMC (β=-0.45, t=-4.04, p=0.0001) as well as between dAIC and PMC (β=-0.39, t=-2.97, p=0.005) ( Figure 5B). Similarly, temporal windows during sleep with greater sleep-task topographic similarity were associated with greater sleep anticorrelation between dPPC and PMC (β=-0.23, t=-1.79, p=0.08) as well as between dAIC and PMC (β=-0.50, t=-3.81, p=0.004). These findings suggest that externally-oriented task-like network states emerge spontaneously during rest and sleep when inter-network anticorrelation is increased.  Figure 4. Similar iEEG infraslow HFB network topography across task, rest and sleep states. A) Locations of seed regions (green for dorsal posterior parietal cortex, blue for posteromedial cortex, red for dorsal anterior insular cortex) and other implanted electrode locations (black) in three subjects included in cross-state wider network-level analysis. B) Spatial correlations in 3 subjects of infraslow HFB functional connectivity values (Fisher-transformed) for all pairs of regions in task versus rest. C) Same as B) but for task versus sleep. D) Same as B) but for rest versus sleep. Green, blue and red data points, respectively, indicate paired targets with the dPPC, PMC and dAIC. Black-outlined yellow, red and orange data points indicate values for dPPC-PMC, dPPC-dAIC, and dAIC-PMC pairs.

Figure 5. Topographic rest-task and sleep-task similarity covaries with the strength of inter-network anticorrelation.
A) Example (from subject S2) channel by channel matrices of <0.1 Hz HFB correlations in 100-second resting state temporal windows with strong (left) and weak (right) dPPC-PMC anticorrelation. Compared to the window with weaker anticorrelation, the window with stronger anticorrelation shows relatively greater topographic similarity with the task template (i.e., the average correlation matrix derived from continuous task performance). B) Across independent 100-second resting state windows, a negative correlation is found between rest inter-network (dPPC-PMC and dAIC-PMC) <0.1 Hz HFB anticorrelation and rest-task topographic similarity (values normalized within subjects and concatenated across subjects). C) Same as B) but for sleep instead of rest.

DAN and SN activations precede DMN deactivations during task performance
Our findings so far, inspired by previous fMRI studies, have considered zero-lag interactions of slow (filtered) activity fluctuations between regions. However, on a finer time scale that is measurable with iEEG, it is possible that inter-network anticorrelated activity is better explained by time-lagged interactions (Raccah et al., 2018). We therefore next sought to determine whether there was a systematic temporal order of iEEG task-evoked HFB responses across networks and whether lagged interactions could be relevant to fluctuations in GradCPT behavioral performance.
Focusing first on participants with simultaneous dPPC-PMC coverage (n=6), we found that the HFB activations in dPPC were significantly earlier than HFB deactivations in the PMC ( Figure   6A). These responses were seen at the individual level during correct omissions but were largely absent during correct commissions ( Figure 6B). In all six participants with relevant coverage, the HFB time-to-peak (TTP) after trial onset for maximum amplitude in the dPPC (M±SD = 772±208 ms) was consistently earlier than that for the peak deactivation in the PMC (1106±222 ms) ( Figure 3C), and this timing difference was statistically significant (p=1.5x10 -6 , Wilcoxon signed rank test on binned subsets of trials). In three of these participants with additional electrode coverage in visual cortical areas, increased HFB power was earlier than those found for the dPPC and PMC ( Figure S3B). Thus, the cross-regional temporal profile was consistent with an expected pattern of information transfer from unimodal to transmodal cortex, with the additional novel finding of a temporal hierarchy within the transmodal cortices themselves (i.e., the DMN deactivations hallmarking the latest stage of processing).
Further analysis of trial-by-trial cross-correlations of inter-regional HFB time series during correct omission trials (i.e., withholding button presses to infrequent target mountain trials) revealed dPPC-PMC and dAIC-PMC anticorrelations that were non-zero-lag i.e., the anticorrelation between the two structures was shifted in time with a significant lead time from dPPC and dAIC to PMC ( Figure S4A). A notable exception was seen in patient S3 (possible explanation provided when we discuss the behavioral significance of shifted anticorrelations). Interestingly, non-zerolag anticorrelations were also seen, to lesser degree, for correct commission trials -i.e., frequent button presses during non-target city trials (Figure S4A), suggesting that temporal coordination between networks was not purely a product of responses to infrequent target stimuli.

Behavioral significance of shifted anticorrelations between DMN and DAN
Previous fMRI evidence has suggested that time windows of greater zero-lag DMN-DAN anticorrelation are associated with better behavioral performance within and across individuals , Thompson et al., 2013, Wang et al., 2016. As our timeresolved iEEG analyses had revealed a shifted (i.e., non-zero-lag) inter-network anticorrelation, we aimed to determine whether these activity lags were similarly behaviorally significant. We

B
hypothesized that time-shifted, but not zero-lag, dPPC-PMC anticorrelation would reflect a subject's level of overall sustained attention that varied across sessions.
To assess the relationship with sustained attention performance (d′), we estimated dPPC-PMC and dAIC-PMC functional connectivity across each task performance session using HFB envelope fluctuations. We expected that behaviorally significant time-lagged interactions would be most apparent in unfiltered, or minimally filtered, but not infraslow signals. We therefore repeated our analyses with unfiltered, 0.1-1 Hz filtered, and <0.1 Hz filtered HFB signals. The 0.1-1 Hz filter range was determined on the basis of our own and others' prior work linking BOLD functional connectivity with HFB power amplitude in this frequency range (Nir et al., 2008, Keller et al., 2013, Foster et al., 2015, Kucyi et al., 2018a. We computed both zero-lag correlations and lag-minimum correlations, defined as the maximum anticorrelation between regions for time series that could be shifted from -2 to +2 seconds ( Figure 8A).
We found that d′ was significantly associated with dPPC-PMC lag-minimum correlation (β 1 =-0.56, t=-3.78, p=0.0006), but not with zero-lag correlation (β 2 =0. 18, t=1.25, p=0.22), across all subjects and sessions for the 0.1-1 Hz HFB range. Specifically, greater dPPC-PMC lag-minimum correlation was associated with better sustained attention (higher d′) (Figure 8B), whereas zerolag anticorrelation was not ( Figure 8C). The direction of this effect was consistent within each patient, with a weakest effect for patient S3. The relationship between d′ and dPPC-PMC lagminimum correlation also remained significant even when using the unfiltered, rather than 0.1-1 Hz filtered, HFB envelope (β 1 =-0.53, t=-2.44, p=0.02), suggesting consistency across multiple time scales of HFB activity. However, as expected, the association was not significant for the infraslow HFB envelope (β 1 =-0.29, t=-0.47, p=0.64; β 2 =0.32, t=0.52, p=0.61). As the distinct response profile of the dPPC in patient S3 compared to the others ( Figure 6A) may have suggested that the recording in this subject could be from a parietal neuronal population that is functionally distinct, we repeated the analyses with S3 omitted. After omitting S3, the strength of the behavioral performance relationship with the dPPC-PMC lag-minimum correlation was similar (β 1 =-0.59, t=-3.71, p=0.0009), and the relationship remained significant when using the unfiltered HFB envelope (β 1 =-0.57, t=-2.64, p=0.01).

Figure 8. Association between time-lagged inter-network anticorrelated activity and behavioral performance across sessions. A) Illustration of how functional connectivity was calculated from continuous HFB 0.1-1 Hz filtered time series. Using an example 20-sec time series (top), a lagged cross-correlation was performed (shifting dPPC relative to PMC and vice versa). (Bottom) The zero-lag correlation was taken as the value with no time shift, whereas the lag-minimum correlation was taken as the minimum value across time shifts. In the main analysis, these metrics were calculated based on whole sessions (typically 6 mins long). B) Within
These findings suggest that attentional state dynamics are reflected in temporally-lagged internetwork (especially DAN-DMN) anticorrelated activity. Our results raise important issues regarding interpretation of the concept of intrinsic inter-network antagonism during active-task and task-free states.

Attentional states and time-resolved anticorrelated inter-network activity
Anticorrelated task-evoked DMN-DAN/SN activity has been found in functional neuroimaging studies across an extensive variety of task conditions involving both stimulus-driven and goaloriented attention (Shulman et al., 1997, Raichle, 2015. States of lapsing attention and mindwandering, which are largely incompatible with sustained external task-oriented attention, have been associated with increased DMN activation, or a lack of DMN suppression (Weissman et al., 2006, Christoff et al., 2009, Kucyi et al., 2016. It has thus been proposed that anticorrelated networks may continuously compete with one another for control of shared computational resources (Sonuga-Barke and Castellanos, 2007. Building on this framework, our finding that task-evoked DAN and SN activations precede DMN deactivations may point toward a causal chain of events that is required for successful deployment of stimulus-driven and goal-oriented attention. During a baseline (or low cognitive demand) state, the DMN maintains control over computational resources (e.g. for imagery associated with internally-oriented cognition). When behaviorally-relevant sensory information is successfully transferred to DAN/SN regions, those regions gain control over resources that the DMN previously had access to, and subsequently the DMN is actively suppressed. An important caveat is that our finding of a systematic temporal order of neuronal population responses does not necessarily imply a causal interaction between DAN/SN activation and DMN suppression, and further study on directional relationships and including other network nodes will be needed. Our results, however, do suggest that the relatively late suppression of the DMN relative to activation in other association networks suggests a temporal hierarchy that may accord well with findings that situate DMN regions as those with longest connectivity paths and furthest geodesic distance from primary sensory regions (Margulies et al., 2016).
The dynamics of anticorrelated activity in relation to task performance have previously been studied largely with group-level fMRI, focused on zero-lag interactions of slow hemodynamic signals. States of greater DMN-DAN/SN anticorrelation have been associated with greater vigilance and behavioral stability (Thompson et al., 2013, Wang et al., 2016, Kucyi et al., 2017. Across individuals, greater baseline BOLD anticorrelation has been associated with lesser response time variability , fluid intelligence , and greater working memory capacity (Keller et al., 2015)-all of which are behavioral measures that may rely on sustained attention. Moreover, attenuated anticorrelation has been found in clinical conditions involving attentional dysfunction , Sripada et al., 2014 as well as in cognitive decline with aging (Keller et al., 2015, Spreng et al., 2016.
Our findings extend those results in multiple ways: First, we provide critical neurophysiological validation for slow DMN-DAN/SN anticorrelated activity during task-related sustained attention at the level of functionally localized neuronal populations and reliably so at the individual level.
Second, we show that electrophysiological anticorrelated activity is associated with momentary changes in sustained attention within individuals. Third, and most importantly, we show that behaviorally-relevant anticorrelated activity involves inter-network lags of up to hundreds of milliseconds that are too short for fMRI to detect. Interestingly, the lagged inter-network anticorrelation magnitude was largely independent of zero-lag infraslow anticorrelation, a finding that accords well with the notion that infraslow brain activity has unique spatiotemporal dynamics compared to faster activity (Mitra et al., 2018). Given that such faster and slower time scales of activity were dissociable, our iEEG analyses may have been sensitive to behaviorally-relevant, time-resolved interactions that would not be detectable with current human neuroimaging methods. Despite the low temporal resolution of functional neuroimaging, inter-network directional and lagged interactions have long been of interest in task and resting states (Nyberg et al., 1996). Due to regional heterogeneity in blood flow dynamics (David et al., 2008), it remains an open question whether advances with accelerated neuroimaging will allow detection of temporallyordered, anticorrelated DMN-DAN/SN activity and its variation over time.

Electrophysiological interactions between DAN, DMN and SN
In the functional neuroimaging literature, the antagonistic relationship between the DMN and DAN -potentially highlighting a competition between internally-oriented and externally-oriented attention -has received intense focus and scrutiny. Though typically lesser emphasized, the SN, also shows negatively correlated BOLD activity with the DMN (as also found within individual subjects here) (Fox et al., 2005, Uddin et al., 2009, Kucyi et al., 2012. Functional neuroimaging evidence indicates that the SN and DAN have dissociable roles in externally-oriented attention. The DAN shows domain-general activation (in tandem with DMN deactivation) during various conditions involving goal-oriented attention (Corbetta and Shulman, 2002). The SN shows activation during detection of salient external stimuli (Downar et al., 2000) and during behavioral errors in fMRI (Neta et al., 2015) and single-unit recordings from key SN nodes (Fu et al., 2018) .
It has been proposed that the SN, and the dAIC in particular, causally facilitates switching between other networks (including DMN and DAN) to reorient attention during salient event detection Uddin, 2010, Uddin, 2015). In partial agreement, a recent application of dynamic causal modeling to resting state fMRI data suggested that the SN and DAN exert intrinsic inhibitory influences on the DMN (Zhou et al., 2018).
Our iEEG results extend these frameworks and confirm the presence of electrophysiological

Intrinsic inter-network anticorrelation
The fMRI-based finding of anticorrelated networks during wakeful rest has led to the notion that there is an intrinsic, state-independent, antagonistic relationship between the DMN and other networks (Fox et al., 2005). Under this framework, the brain may continuously shift between states that draw, respectively, from internally-and externally-oriented sources of information (Buckner et al., 2008, Honey et al., 2018. However, the concept has remained controversial, in large part due to technical limitations of fMRI (Murphy and Fox, 2017). Infraslow resting state BOLD anticorrelations become introduced into data following preprocessing with global signal regression (Murphy et al., 2009), but anticorrelations have also been detected in the absence of global signal regression and with alternative noise-correction strategies (also found within our cohort) (Fox et al., 2009, Chai et al., 2012. Though anticorrelation of infraslow inter-network activity has been recovered in computational models (Deco et al., 2009) We investigated the electrophysiology of inter-network anticorrelation with several advances: 1) We studied functionally localized DMN-DAN and DMN-SN region pairs that were identified based on task-evoked responses; 2) We studied extended (i.e., 24-48 minutes) recordings across task, rest and sleep states; and 3) We investigated infraslow HFB activity, which closely resembles the typical temporal scale of BOLD activity studied during wakeful rest. Based on these approaches, we reliably observed anticorrelated activity during continuous task performance but not during Importantly, our findings do not conclusively rule out the possibility that intrinsic antagonistic relationships exist in the brain. We focused here on network nodes where electrode coverage was available, often only within a few sites per subject (especially in depth electrode cases). Consistent with previous studies (He et al., 2008, Keller et al., 2013, Ramot et al., 2013, Foster et al., 2015, Hacker et al., 2017, Kucyi et al., 2018a, we found that the spatial topography of functional connectivity patterns remained similar across task, rest and sleep states and between iEEG and within-individual resting state fMRI. We also found that even when dPPC-PMC and dAIC-PMC anticorrelations were not found, inter-network correlations between these pairs were typically within the lowest percentiles compared with all other region pairs. A possible explanation is that DMN-DAN/SN pairs have an increased propensity to exhibit snippets of spontaneous anticorrelated activity. Indeed, we found that rest and sleep anticorrelations were associated with the expression of 'task-like' network states. It is thus possible that variations in inter-network activity coordination could be dependent on cognitive state. Though some form of behavioral validation would be needed to confirm, our task findings of increased anticorrelation with greater sustained attention support this possibility. It should be acknowledged, however, that variations in functional connectivity can arise due to sampling variability . There remains an ongoing debate over whether spontaneous temporal fluctuations in functional connectivity within and between networks are behaviorally significant (Liegeois et al., 2017, Kucyi et al., 2018b.

Conclusion
Our findings establish a behavioral significance of systematic temporal lags underlying anticorrelated inter-network activity. Additionally, our work suggests that if intrinsic antagonistic inter-network relationships exist, their expression may vary across time and as a function of behavioral state. This knowledge is critical for the interpretation of task and resting state functional neuroimaging studies and for understanding the basis of changes in inter-network relationships in health, aging and disease. The temporally ordered inter-network interactions identified here point toward the possible capacity for causal influences, a topic that requires further study with neuromodulatory techniques such as direct brain stimulation.

Subjects
Data from seven human subjects (S1-S7) who were undergoing neurosurgical treatment for refractory focal epilepsy were included in analyses reported here. Data from S1-S6 were collected at Stanford University Medical Center, whereas data from S7 were collected at Beijing Tian Tan Hospital (age range: 19-34, 4 females and 3 male, all right-handed; see Table S1 for full demographic and other details). Subjects were implanted with intracranial electrodes (Adtech Medical Instruments) that either were depth electrodes placed stereotactically within both hemispheres (S2, S3, S7) or one hemisphere (S4), subdural electrodes arranged in grid and strip configurations over one hemisphere (S5, S6), or a mixture of both (S1, who had all depth electrodes except for one strip covering ventral temporal cortex  Table S1 for number of electrode contacts within each of these regions for each patient). The seven patients included here were selected from a cohort of 32 patients (10 at Stanford, 22 at Beijing) who participated in the cognitive task procedures described herein. Out of those patients, 10 had simultaneous coverage in the regions of interest. Out of those 10, 3 were excluded for the following reasons: 1) irregular signals in the PMC (according to criteria defined under Intracranial EEG: Data Preprocessing); 2) encephalomalacia found in the occipital lobe; and 3) failure to record signals from a region of interest due to hardware failure.

Intracranial EEG Data Acquisition
Intracranial EEG (iEEG) recordings were performed at bedside of the subject's private clinical suite. For Stanford patients (S1-6), data were recorded either via a multichannel research system (Tucker Davis Technologies, Alachua FL, USA) or a clinical monitoring system (Nihon Kohden, Tokyo, Japan). For the Beijing patient (S7), all data were recorded with a Nihon Kohden system.
The sampling rate was either 1000 or 2000 Hz for all task sessions and for some wakeful rest sessions, whereas the sampling rate was either or 500 or 1000 Hz for wakeful rest and sleep sessions (see Table S2 and contacts were circle-shaped with diameter of 2.3 mm in the exposed area of recording and interelectrode spacing of 5-10 mm. For the Beijing patient, depth electrode contacts (HKHS Healthcare, Beijing, China) had a contact length of 2 mm, diameter of 0.8 mm, and inter-electrode spacing of 1.5 mm. During recording, the iEEG signals were referenced to the most electrographically silent channel outside of the seizure focus. The total number of electrode sites ranged from 60 to 210 (Table S1).

Continuous Performance Task Sessions
The gradual-onset continuous performance task (GradCPT) (Esterman et al., 2013) was administered in multiple (range: 4 to 8) sessions in each patient, each lasting 4-8 minutes ( Table   S2). The number of sessions obtained within each patient depended on time available for research testing in the clinical environment, which varied across patients. The task was administered at bedside via a laptop (running Windows 10 Pro and Windows 8.1, respectively, in Stanford and Beijing) with its screen positioned ~70 cm from the patients' eyes at chest level.
Stimuli were presented using Psychophysics Toolbox (Brainard, 1997) in Matlab R2016b (MathWorks, Natick MA, USA). An RTBox device (Li et al., 2010) was used to send transistortransistor logic pulses to an empty channel on the EEG montage to mark the onset times of each stimulus.
During task performance, grayscale visual images of either city or mountain scenes appeared within round frames (with white background) and gradually transitioned from one to another for the duration of the task. Each transition lasted 800 ms. Using linear pixel-by-pixel interpolation within each trial, image coherence began to gradually increase from time zero (minimum coherence) until 400 ms (maximum coherence) before gradually decreasing back to minimum coherence (at 800 ms). Scenes were presented randomly with 10% mountain and 90% city, but the same scene could not repeat on consecutive trials. In one subject (S7), two additional runs were performed with 25% mountain and 75% city rates, and these runs were included in functional localization analyses but were excluded from analyses of task performance versus neural activity due to potential variation in task difficulty. Subjects were instructed to press the space bar on the laptop upon noticing each city appearing but to withhold response when noticing a mountain appearing. Subjects were asked to perform their best and to keep going when they noticed themselves making an error. Each task session began with a 20 second baseline period in which the patient was instructed to fixate on a blurred mask stimulus (same size as scene stimuli) and get ready to begin. Subjects performed with their dominant hand, except for S4 who performed with their non-dominant hand due to discomfort of the dominant hand.

Rest and Sleep Sessions
Resting state and sleep sessions were obtained and analyzed within three patients (S1, S2, and S6). Resting state recording sessions were obtained either with instructions (using a research system) or without instructions (using the clinical monitoring system) (see Table S3 for summary).
In instructed rest sessions, subjects were asked to relax and not think of anything in particular while either keeping eyes closed or open, including select sessions where subjects fixated on a central fixation cross on the laptop screen. To increase the total duration of resting state recordings, we reviewed video EEG during daytime hours and clipped periods (each at least several minutes long) when the subject appeared relaxed, awake and not engaged in conversation, clinical procedures, or with television/handheld electronic devices.
To obtain sleep recordings, we reviewed video EEG of recordings during nocturnal hours and clipped periods (each at least several minutes long) when the subject appeared to be asleep. A total duration was obtained that was similar to that obtained during task and resting states (see Table S4 for summary). Though precise sleep staging data were not available, we inspected lowfrequency power at iEEG channels of interest during wakeful and sleep recordings, as done previously (Foster et al., 2015). This allowed us to identify sleep-related features such as increased delta power, attenuation of alpha/theta oscillations and presence of sleep spindle-like oscillations at ~15 Hz ( Figure S5).

Anatomical Localization of Electrode Contacts
We used the iElvis pipeline (Groppe et al., 2017) for anatomical localization of electrode contacts.
First, we processed and reconstructed the T1 scan using Freesurfer v6.0.0 (recon-all command) (Fischl et al., 1999). We then aligned the post-implant CT image to the pre-implant T1 scan using a rigid transformation (6 degrees-of-freedom, affine mapping), and we inspected the quality of the registration. Using BioImage Suite (Papademetris et al., 2006), we manually labeled each electrode contact location on the T1-registered CT image. For subdural electrode cases (S5, S6), we then projected the electrode coordinates to the leptomeningeal surface and applied correction for post-implant brain shift, using previously described methods (Dykstra et al., 2012). In stereotactical EEG cases, minimal post-implant brain shift is expected, and thus no further adjustment was made. The electrode coordinates obtained from these approaches were used for visualization and for fMRI analyses.

Intracranial EEG: Data Preprocessing
Data from each iEEG recording session were preprocessed similarly for task, rest and sleep sessions using an identical pipeline to that described in our previous work (Kucyi et al., 2018a).
The procedures drew from tools in the Matlab-based LBCN preprocessing pipeline (https://github.com/LBCN-Stanford/Preprocessing_pipeline), SPM12 (Kiebel and Friston, 2004), and Fieldtrip (Oostenveld et al., 2011). For task runs, the recording was first cropped to retain data only within the baseline and task performance periods. Notch filtering was performed to attenuate power-line noise (band-stop between 57-63, 117-123, and 177-183 Hz for data from Stanford, and band-stop between 47-53, 97-103, and 147-153 Hz) for data from China. We then re-referenced the signal from each channel to the common average signal across all channels, with the following channel types excluded from the common average those that: a) showed pathological activity during clinical monitoring (as noted by a neurologist); b) were manually labeled as clear outliers on power spectra plots of all channels; c) had a variance greater or lesser than five times the median variance across all channels; or d) had greater than three times the median number of spikes across all channels, with spikes defined as 100 µV changes between successive samples. We then performed time-frequency decomposition using a Morlet wavelet transform with frequencies of interest log-spaced between 1 and 170 Hz (38 total values).
To normalize the distributions of power amplitude estimates, for each frequency of interest, we rescaled each time sample by the log ratio of the whole session's power amplitude time series.
This rescaling step accounted for the band-specific 1/f decline of the power spectrum (Miller et al., 2007). Subsequently, we performed averaging of power amplitude estimates within seven frequency bands, including δ (1-3 Hz), θ (4-7 Hz), α (8-12 Hz), β1 (13-29 Hz), β2 (30-39 Hz), γ (40-70 Hz) and HFB (70-170 Hz). We then visually inspected the HFB time series at each channel in each session, and we excluded channels that showed irregular, spikey or pathological activity (that may have been otherwise missed in our inspection/exclusion prior to time-frequency decomposition). Additionally, we manually reviewed the anatomical locations of electrodes and removed from analysis those channels that were located outside of the brain or largely within white matter. For each participant, we computed the average number of channels retained across sessions.

Resting fMRI Preprocessing
All fMRI volumes were manually inspected for possible artifacts and head movements. The mean relative head displacement values for subjects (S1-S5) were 0.05, 0.03, 0.05, 0.03, and 0.09 mm.
Subsequently, the segments were eroded to retain the top 198 cm 3 and top 20 cm 3 of voxels with highest probability of being WM and CSF, respectively (Chai et al., 2012). This erosion was performed to minimize contamination of GM signal within the WM and CSF volumes. We then regressed out of each voxel the following: mean global signal, mean WM signal, mean CSF signal, and 6 motion parameters obtained with MCFLIRT. Finally, spatial smoothing (6mm full width at half maximum kernel) and temporal filtering (0.01-0.1 Hz) were performed. Importantly the above-described pipeline involved global signal regression, a procedure that remains controversial, especially in the context of anticorrelated brain networks (Murphy et al., 2009, Murphy andFox, 2017). We therefore also present results from an alternative preprocessing pipeline including ICA-AROMA (Pruim et al., 2015) and without global signal regression. For this pipeline, we performed deletion of the first 4 volumes, brain extraction, motion, spatial smoothing (6mm full width at half maximum kernel), and nonlinear registration across fMRI, T1, and MNI152 standard spaces. We then applied independent components analysis (ICA) with FSL's MELODIC and automatic dimensionality estimation, and components that were classified as noise were regressed out of each voxel, as described previously (Pruim et al., 2015). Bandpass temporal filtering (0.01-0.1 Hz) was then applied.

Anatomical Classification of Electrode Contacts
We conducted region-of-interest (ROI) based analyses of the dorsal parietal cortex (dPPC), posteromedial cortex (PMC) and dorsal anterior insula (dAIC) due to their well-described memberships within the dorsal attention, default mode and salience networks, respectively. While other nodes of these networks were of interest, we focused on these three nodes in part due to practical considerations, as electrode coverage was based solely on clinical decision-making. We classified channels as being within the dPPC, PMC or dAIC based on individual-level anatomy reviewed on 3D T1 volumes and cortical surface reconstructions (see Table S1 for summary of number of channels per ROI identified in each patient).
As both the superior parietal lobule and intraparietal sulcus (IPS) have been linked to the dorsal attention network Shulman, 2002, Fox et al., 2006), we considered these adjacent areas within the parietal lobe as a single ROI, which we term dPPC. We defined the PMC in a similar fashion as in previous work (Parvizi et al., 2006, Dastjerdi et al., 2011, Foster et al., 2015. The PMC included areas posterior to the post-central sulcus within the posterior cingulate cortex (within BA 23a and 23b), retrosplenial cortex (BA 29/30), and medial parietal cortex/precuneus (BA 31 and 7m). These areas were bounded by the marginal branch of the cingulate sulcus (dorsally/anteriorly) and by the parieto-occipital sulcus (posteriorly).
The dAIC was defined based on boundaries and landmarks defined previously (Ture et al., 1999, Naidich et al., 2004 to demarcate areas that corresponded largely to the agranular anterior insular zone (Mesulam and Mufson, 1982). This included the accessory gyrus of the insula, and/or portions of the anterior, middle and posterior short gyri of the insula that were superior to the inferior-most point of the short insular sulcus. These dorsal anterior subregions of the insula have been consistently linked with the salience, or cingulo-opercular, network (Seeley et al., 2007, Kelly et al., 2012.
For a supplemental analysis, we also defined higher-order visual cortical areas, including the lingual gyrus, fusiform gyrus and inferior temporal gyrus.

Functional Localization of Task-responsive iEEG Channels
After identifying all electrode contacts within each ROI, we defined functionally responsive channels during GradCPT performance. Specifically, we assessed evoked HFB power during correct omissions (withheld behavioral responses) to rare, target trials (mountain scenes) relative to correct commissions (behavioral responses) to frequent, city trials. We also compared HFB responses during correct omission versus commission error (incorrect behavioral response) trials.
Based on replicated findings from previous fMRI studies (Esterman et al., 2013, Fortenbaugh et al., 2018, and on the known association between BOLD activity and electrophysiological HFB activity (Logothetis et al., 2001, Mukamel et al., 2005, Nir et al., 2007, Hermes et al., 2012, we expected that HFB power would show an increase at dPPC and dAIC as well as a decrease at PMC during correct omissions and commission errors. For this analysis, we minimally smoothed the HFB power amplitude time course within each session using a 50-ms Gaussian window. We then extracted the HFB time series from windows surrounding each mountain trial, with each window starting at 800 ms prior to mountain scene onset (start of fade-in) and ending at 1600 ms after the onset. We excluded from this analysis mountain trials that were preceded by other mountain trials. We also extracted time windows with the same boundaries around correct commission responses (city trials with button presses). For these correct commission trials, we extracted only those that were both preceded and followed by other city trials (and both with correct responses) to avoid potential contamination with responses evoked by rare mountain scenes. Among these retained correct commission trials, we deleted a random subset of trials within each session such that the remaining subset included a total number of trials that matched the number of correct omission trials within the same session. All GradCPT sessions were included in this analysis except for two runs in one patient (S1), where excessively poor performance was found (>75% commission error rate), suggesting that the patient may not have been attending to target events.
To assess significance of HFB responses during correct omission compared to correct commission trials as well as correct omission compared to commission error trials, we adopted a nonparametric cluster-based permutation test as implemented in Fieldtrip (Maris and Oostenveld, 2007)  showed the cluster with the largest effect among all significant clusters (i.e., lowest significance probability value), and the peak-responsive channel within each ROI was selected for focused analyses described below. To evaluate the specificity of effects to HFB relative to other frequency ranges, we inspected spectrograms of mean responses during correct omission trials based on rescaled power amplitude estimates in the range of 1-170 Hz (Figure 1) (for visual purposes, smoothing with a Gaussian kernel was applied).

Multiple Kernel Learning Analysis of iEEG Responses
To comprehensively assess the possible contributions of different frequency bands of activity to task-evoked iEEG responses, we performed a multiple kernel learning (MKL)-based analysis. The MKL approach is a machine learning method that can be applied to classifying iEEG task conditions by including multiple frequency bands of activity as well as multiple channels in a single model (Schrouff et al., 2016). We used MKL in the PRoNTo toolbox (Schrouff et al., 2013) to classify correct omission versus correct commission trials.
For each subject, all channels that were anatomically identified as being within one of the ROIs (PMC, dPPC, dAIC) were included in the model. In a full model, channels from all ROIs within a given subject were included, and in additional single ROI models, channels only from given ROIs were included. For each trial, power amplitudes from each channel were extracted between time 0 to 1500 ms after trial onset for seven frequency bands (δ,θ,α,β1,β2,γ,HFB). The number of trials was matched between conditions (as described above). Model features were defined as kernels, or pair-wise similarity matrices across the time series of all trials, which were constructed for each channel and frequency band (i.e., the number of kernels per subjects was m x 7, with m being the number of channels). Each kernel was normalized and mean-centered to ensure that modeling was not influenced by the scale of each kernel.
We then applied MKL, using a support vector machine to define a decision boundary to discriminate between correct omission and correct commission trials. As in previous work (Schrouff et al., 2016), model parameters were optimized to determine the decision boundary for each kernel, and decision boundaries were weighted by a parameter d m to define a global decision boundary. We used a 10-fold cross-validation scheme: in each fold, training was performed 90% of trials, and testing was performed on the 10% left out trials (with a different 10% left out on each fold). Model accuracy was obtained as the average balanced accuracy (average of class accuracies) across folds. During cross-validation, the soft-margin parameter, C, was optimized by considering values 0.01, 0.1, 1, 10, 100 and 1000. A nested cross-validation was performed where the value of C leading to highest model performance in the inner crossvalidation was selected, and that C value was used to estimate performance in the outer crossvalidation. Statistical significance of model accuracy within each subject was assessed using 1000 permutations of the training labels to generate an accuracy distribution, and p values less than 0.05 for the true value were considered as significant. To evaluate the contributions of different frequency bands of activity to model performance, for each fold we calculated the sum of d m values across channels for each of the seven frequency bands. We then calculated the mean of those sums across the ten folds.

Resting fMRI: Seed-based Functional Connectivity
We transformed the within-subject coordinates obtained from electrode localization to fMRI space (using the previously computed linear transform). We then extracted the BOLD time series from seed regions defined as 6-mm radius spheres at electrode locations of interest (i.e., peakresponsive channels, defined based on criteria described above). Using a general linear model (implemented in FSL) for each seed region separately, the demeaned the BOLD time series was entered as a regressor. We then projected the resulting volume-map z-scores obtained at each voxel to vertices on the cortical surface in Freesufer. To provide a comparison between individual-level networks found with this approach and the standard networks found in healthy populations, we registered the DMN, DAN and SN templates from the 7-network Yeo parcellation (Yeo et al., 2011) from standard (fsaverage6) to individual surfaces.
In the two patients who did not undergo fMRI, we registered the Yeo parcellation to individual cortical surfaces to determine whether peak-responsive channels were within the DMN, DAN or SN. As one of these two patients had depth electrodes (rather than subdural surface recordings), raw electrode coordinates were in volume space rather than in Freesurfer cortical vertices. Thus we snapped the coordinates to the nearest cortical vertex and determined its Yeo network identity.

Intracranial EEG: Comparison of Task, Rest and Sleep States
We compared HFB functional connectivity across states (task, rest and sleep) for dPPC vs. PMC and dAIC vs. PMC, using the functionally localized peak-responsive channels for each subject within ROIs. To increase the statistical power of this analysis, we divided each task, rest and sleep recording session into 100-second, non-overlapping windows (deleting the remaining last part of each recording). The 100-second window length was selected so that a slowest frequency component of 0.01 Hz, which is commonly assumed to be relevant to infraslow functional connectivity (Leonardi and Van De Ville, 2015), was retained. Within each window, we computed functional connectivity as the Fisher-transformed inter-regional correlation of the infraslow (<0.1 Hz) HFB envelope. Using these functional connectivity values, we performed a three-way analysis of variance (ANOVA) within each subject with factors Task, Rest and Sleep (significance set at two-tailed p<0.05). We performed post-hoc two-sample t-tests to compare functional connectivity between pairs of states and one-sample t-tests to assess whether non-zero functional connectivity was reliably detected within each state (significance set at two-tailed p<0.05). We also repeated the same analyses using power amplitude correlations in six other frequency bands other than HFB (δ, θ, α, β1, β2, and γ).
We additionally performed a wider network-level functional connectivity analysis to compare the spatial topography of HFB inter-regional correlations across states. For this analysis, all channels within each subject that were retained after preprocessing were included. Using the peakresponsive channels within the dPPC, PMC and dAIC as seed regions, we computed the interchannel time series correlations (Fisher-transformed) within each session and then averaged these values across sessions for each state (i.e., generating a vector for each seed region to all target regions). To avoid potential spurious correlations due to volume conduction, we excluded target channels that were immediately neighboring the seed channel (on the same depth probe or subdural strip). We then concatenated vectors across seed regions within subjects and deleted values from redundant pairs. We then correlated the obtained vectors from the three states (task, rest sleep) with one another as an index of the inter-state spatial similarity of HFB functional connectivity (significance set at p<0.05). We also repeated these analyses of wider network-level functional connectivity using power amplitude correlations in six other frequency bands other than HFB (δ,θ,α,β1,β2,and γ).
We performed additional analyses to test the hypothesis that task-like topographic network patterns would be found during temporal windows in which rest and sleep infraslow HFB anticorrelations emerged. For this analysis, we computed a 'task template' topographic pattern, defined as the mean correlation matrix (Fisher-transformed and using all channels retained for analysis) across all independent 100-second windows during task performance. For each rest and sleep 100-second window, we computed the correlation matrix (Fisher-transformed) between the same channels. We then correlated each rest and sleep matrix with the task template. For the two ROI pairs of interest (dPCC-PMC and dAIC-PMC) and the two behavioral states (rest and sleep), we then performed linear mixed model analyses, as implemented in the R environment (Baayen et al., 2008). Across all subjects and temporal windows, subject was entered as a random effect, similarity to task template was entered as dependent variable, and ROI pair correlation was entered as fixed effect. For these analyses, task template similarity and ROI pair correlation values were within-subject normalized (by mean and standard deviation). Significance was set at p<0.05 (Sattherthwaite's approximation) (Luke, 2017).

Intracranial EEG versus BOLD Functional Connectivity Analysis
For each state recorded in iEEG data (task, rest, sleep) we compared the spatial topography of HFB functional connectivity with that of BOLD functional connectivity within subjects. For this analysis, we extracted the BOLD time series from 6-mm radius spheres surrounding each channel's coordinates in fMRI space (as described above). When this procedure resulted in overlapping voxels between regions (e.g. neighboring channels), we deleted those voxels from each region before extracting the time series. Using the locations of the peak-responsive channels within the dPPC, PMC and dAIC as seed regions, we subsequently computed the interregional BOLD correlation (Fisher-transformed) with all target regions (excluding those regions located at channels that were excluded from iEEG analysis). We then concatenated vectors across seed regions within subjects and deleted values from redundant pairs. We then correlated the BOLD vectors with the iEEG vectors from each state (task, rest, and sleep) with regions/channels aligned across modalities (significance set at p<0.05).

Intracranial EEG: Time to Peak Estimation
We performed time-to-peak (TTP) estimation of HFB responses during correct omission trials on peak-responsive channels within the dPPC, PMC and dAIC identified within each subject. Within a time window ranging from +200 to +1500 ms after trial onset (i.e., beginning of mountain scene fade-in), we identified the maximum peak time point for dPPC and dAIC and the minimum peak time point for PMC based on the average HFB response across trials ( Figure 3C). However, because only 6, and 4, subjects were included in these analyses for dPPC vs. PMC and dAIC vs.
PMC, respectively, we performed additional data processing so that more powerful statistics could be performed. Specifically, we divided all trials into equally-sized bins of trials, selecting 8 total bins per subject (i.e. the maximum number of sessions per subject). In each bin we computed TTP as described above. We entered TTP within all bins and across all subjects into two-sided Wilcoxon signed rank tests comparing a) dPPC versus PMC; and b) dAIC versus PMC (significance set at p<0.05).
In a complementary approach to TTP analysis, we performed cross-correlations between channels' time series, with shifts ranging from -2 to 2 seconds. These cross-correlations were performed on single correct omission (mountain) trials, and the mean across cross-correlation across trials was plotted ( Figure S4). We repeated this cross-correlation analysis for all correct commission (city) trials that were preceded and followed by other city trials.

Behavioral Analysis
We used sensitivity (d′) as a measure of task performance, based on signal detection theory (Macmillan and Creelman, 2004), within each session: where Z(p) is the inverse of the cumulative distribution function of the Gaussian distribution. Thus, the higher the d′ value, the higher the overall accuracy of behavioral performance (based on responses to both cities and mountains). Evidence indicates that d′, based on GradCPT performance, is a generalizable measure of sustained attention (Fortenbaugh et al., 2015, Rosenberg et al., 2016. Because of the fast pace of the task and the overlap of stimuli across adjacent trials, key presses were assigned to trials using a previously described iterative algorithm (Esterman et al., 2013, Kucyi et al., 2016. Presses were assigned relative to the beginning of each image transition. For trials in which the reaction time (RT) was highly deviant (before 70% image coherence for current trial, or after 40% coherence for the following trial), the following criteria were used for trial assignment: 1) If the previous or current trial had no response, the press was assigned to the trial in which the response occurred; 2) If both adjacent trials had no response, the press was assigned to the trial closest in time (excluding cases where the trial was a mountain image). 3) If multiple presses could be assigned to a given trial [based on 1) and 2)], the fastest RT was assigned to that trial.

Task-based iEEG Functional Connectivity Analysis
Prior to functional connectivity analysis, we applied a bandpass temporal filter (butterworth, 4 th order) to the unsmoothed HFB envelope, retaining frequencies between 0.1-1 Hz (Nir et al., 2008, Keller et al., 2013, Kucyi et al., 2018a. We performed additional analyses of the HFB envelope, based on no filtering (minimally smoothed, as described above) as well as lowpass-filtering of the time series to the infraslow (<0.1 Hz) range. Task-based functional connectivity analysis was applied for peak-responsive channels in the dPPC, dAIC and PMC, using two metrics: 1) Zero-lag FC: the zero-lag correlation between the channels' time series, and 2) Lag-minimum FC: the minimum correlation (i.e., the greatest anticorrelation) among cross-correlations between the channels' time series, with shifts ranging from -2 to 2 seconds. We then applied a Fisher r-to-z transformation to these values.
To compare functional connectivity with session-to-session variability in behavioral performance (d′) across subjects, we first normalized functional connectivity and d′ values within subjects (i.e., for each session, we subtracted out the mean and then divided by the standard deviation of values across sessions). We then performed linear mixed effect model analyses with subjects entered as random effects, behavioral performance as dependent variable and functional connectivity (zero-lag FC and lag-minimum FC) as fixed effects. Significance was set at p<0.05 (Sattherthwaite's approximation). For visual display purposes, we plotted normalized functional connectivity versus (d′) values, including regression lines illustrating group-and individual-level data (Figures 3E and 3F). Only GradCPT sessions that included 90% city and 10% mountain rate (Table S2) were included in these analyses so that task difficulty was matched across sessions.

Supplemental Figure Legends
Figure S1. Functional localization of dPPC, PMC and dAIC in all subjects (except S2; see Figure  1). (Left) Anatomical locations of electrode contacts implanted in the dPPC, PMC, and dAIC that showed peak-responsive high-frequency broadband         Table 3. Summary of rest sessions conducted.