Disambiguating the role of blood flow and global signal with Partial Information Decomposition

In resting state functional magnetic resonance imaging (rs-fMRI) a common strategy to reduce the impact of physiological noise and other artifacts on the data is to regress out the global signal using global signal regression (GSR). Yet, GSR is one of the most controversial preprocessing techniques for rs-fMRI. It effectively removes non-neuronal artifacts, but at the same time it alters correlational patterns in unpredicted ways. Furthermore the global signal includes neural BOLD signal by construction, and is consequently related to neural and behavioral function. Performing GSR taking into account the underlying physiology (mainly the blood arrival time) has been proved to be beneficial. From these observations we aimed to: 1) characterize the effect of GSR on network-level functional connectivity in a large dataset; 2) assess the complementary role of global signal and vessels; and 3) use the framework of partial information decomposition to further look into the joint dynamics of the global signal and vessels, and their respective influence on the dynamics of cortical areas. We observe that GSR affects intrinsic connectivity networks in the connectome in a non-uniform way. Furthermore, by estimating the predictive information of blood flow and the global signal using partial information decomposition, we observe that both signals are present in different amounts across intrinsic connectivity networks. Simulations showed that differences in blood arrival time can largely explain this phenomenon. With these results we confirm network-specific effects of GSR and the importance of taking blood flow into account for improve denoising methods. Using GSR but not correcting for blood flow might selectively introduce physiological artifacts across intrinsic connectivity networks that distort the functional connectivity estimates.


Introduction
In recent years there has been increasing interest in the use of resting state functional magnetic resonance imaging (rs-fMRI) in neuroimaging research. A popular approach in rs-fMRI is to map the functional architecture of the human brain using patterns in resting-state or intrinsic correlations (Fox and Raichle, 2007). The correlations of low frequency oscillations present in the blood oxygenation level dependent (BOLD) signal reflect the functional connectivity (FC) between different brain regions. Regions that then show a high mutual correlation are referred to as a resting-state network (RSN) or intrinsic connectivity network (ICN). In rs-fMRI studies a common preprocessing step before estimating FC is global signal regression (GSR), where a global time course is regressed out of the data. The global signal (GS) is obtained by averaging the resting-state time courses over the entire brain (Desjardins et al., 2001). The GS is often thought of as a mixture of processes that confound the BOLD fMRI signals (artifacts). Based on the assumption that processes that are globally spread across the brain cannot be linked to neuronal activation, it could be beneficial to remove them to denoise the data. Indeed, fluctuations in the GS have been linked to physiological fluctuations, mainly respiratory effects, head motion, hardware scanner related effects and vascular effects (Murphy and Fox, 2017;Power et al., 2017b).
However, in addition to these non-neuronal confounds, GSR has also been related to neuronal fluctuations. Sch€ olvinck and colleagues showed that the fMRI signal calculated over the entire cerebral cortex in monkeys showed positive correlations with the spontaneous fluctuations in the local field potentials in a single cortical site (Scholvinck et al., 2010). Others have shown a relationship between the global signal and EEG measures of vigilance and broadband electrical activity using simultaneous EEG and fMRI (Wen and Liu, 2016;Wong et al., 2013). These issues are well summarized in (Liu et al., 2017) which set to disambiguate the nuisance and the information component of the GS by looking at its origin and associations with other neuroimaging measures. More recent studies have demonstrated that the activity of the basal forebrain is intimately linked to cortical arousal  and GS (Turchi et al., 2018), showing that inactivation of the basal forebrain leads to increased global spontaneous fluctuations. The GS has been explicitly connected to these fluctuations, either by predicting the nature of the quasiperiodic patterns of large-scale brain activity (Yousefi et al., 2018), or by encoding the transitions between them (Gutierrez- Barragan et al., 2019). Using calcium-based imaging in mice Matsui and colleagues identified global propagating waves of activity in the neocortex of mice, which points to the existence of global neuronal signals (Matsui et al., 2016). In sum, evidence demonstrates that the GS is a mixture of neuronal and non-neuronal components, but it's unclear how much each component contributes to the GS (Uddin, 2017). As a result, GSR may remove fluctuations of neuronal origin, which could induce errors in FC estimates . After all, the whole brain can be thought of as the roughest parcellation in terms of ROIs or Intrinsic Connectivity Networks, and in this sense it contains all neural/behavioral correlates, as well as all the nuisances.
The practice of GSR has been debated since Murphy and colleagues disputed the finding of anticorrelations between the default-mode network (DMN) and the task positive network (TPN) reported by Fox and colleagues (Fox et al., 2005;Murphy et al., 2009). By using a mathematical argument Murphy and colleagues showed that the anticorrelations found were an artifact introduced by GSR, and in the absence of this preprocessing step the DMN and TPN were positively correlated. Mathematically GSR shifts the distribution of FC estimates to be centered around zero, thereby inducing the existence of both positive and negative correlations (Murphy et al., 2009). Similarly, other studies also observed anticorrelations only when GSR was applied (Ibinson et al., 2015;Weissenbacher et al., 2009), while others have found anticorrelations without applying GSR, but after ingestion of caffeine (Wen and Liu, 2016;Wong et al., 2013), physiological noise correction (Chang and Glover, 2009) or component-based noise reduction (Chai et al., 2012), suggesting that anticorrelations might not be an artifact of the GSR preprocessing step. The nature of the anticorrelations as a true reflection of FC versus an artifact of the GSR technique is not clear nor agreed upon. But most evidence seems to point to a reduction in positive correlations and an increase in spurious anticorrelations. Moreover, using a modeling approach it has been shown that GSR does more than "just" inducing anticorrelations, also altering the underlying correlation structure in unpredictable ways. Saad and colleagues showed in a group comparison study that applying GSR alters short-and long-range correlations within a group, leading to spurious group differences in regions that were not modeled to show true FC differences (Saad et al., 2012).
Despite the critiques a consensus on the use of GSR has not been reached, and it remains a popular denoising method in rs-fMRI studies (Murphy and Fox, 2017;Power et al., 2017b). A recent wave of studies further proved its usefulness in reducing the impact of artifacts, particularly those resulting from participant head motion. Studies comparing the effectiveness of different denoising approaches show that only combinations of preprocessing pipelines that include GSR can minimize the effects of motion Lydon-Staley et al., 2019;Parkes et al., 2018;Satterthwaite et al., 2017), temporally lagged artifacts (Byrge and Kennedy, 2018) and respiration (Power et al., 2017a).
Recent work by Erdo gan and colleagues proposed an improvement to the standard GSR method called dynamic GSR (dGSR) which improves the sensitivity of FC estimates (Erdo gan et al., 2016). dGSR effectively deals with systemic low frequency oscillations (sLFO's) that are non-neuronal of origin and oscillate at frequencies of interest in rs-fMRI (~0.01-0.1 Hz). Here we refer to sLFO as fluctuations of clear non-neuronal origin (as defined in the cited studies, see Tong et al., 2019 for a recent comprehensive review), as opposed to quasiperiodic patterns (QPPs) which reflect neural activity. These intrinsic signals appear to have a vascular origin that travel with the blood through the body (Tong et al., 2013;Tong and Frederick, 2012) and brain (Tong and Frederick, 2010). The sLFO's appear in the arteries before propagating through the brain with the cerebral blood flow from large arteries to large veins (Tong et al., 2018). The spatial-temporal pattern of the arrival time of these sLFO's are consistent with blood flow circulation patterns obtained with dynamic susceptibility contrast (DSC) bolus track imaging (Tong et al., 2017). The method is based on taking the temporal information of blood arrival time into account. For each voxel an optimally delayed version of the global signal is regressed out, which effectively removes the systemic noise of these sLFO's. As a result, the strength and specificity of FC estimates is improved compared with standard GSR. In addition, in further research they show that BOLD signal from major arteries (internal carotid artery) and prominent veins (superior sagittal sinus, SSS and internal jugular vein, great vein of Galen) are highly correlated with the GS (Liu et al., 2017;Tong et al., 2018), supporting the importance of a vascular component in the GS. From this group of studies, it's evident that a large portion of the GS has a macro-vascular origin, and that considering the temporal blood flow information improves de-noising methods. Indeed, taking information from the vessels in de-noising approaches is important, as recent work has shown that the dynamic effects of vessel information has complex consequences on BOLD responses (Kay et al., 2019). Some recent efforts have been made to merge perfusion and BOLD to take blood into account in FC studies (Tak et al., 2015;Cohen et al., 2017).
In a recent correspondence, three directions forward are suggested to further evaluate the use of GSR (Power et al., 2017a(Power et al., , 2017bUddin, 2017). One of these directions raises the concern that the field lacks empirical studies that focus on the effect of GSR in high dimensional fMRI data. In this study, we attempt to address this issue by applying GSR to high dimensional empirical fMRI data, using an information theory and correlation approach to investigate how FC estimates are affected by its use. We are interested in a more detailed evaluation of which brain regions and intrinsic connectivity networks in the connectome are affected by GSR. In line with the work from (Gotts et al., 2013;Saad et al., 2012), we expect that the large-scale brain connectivity pattern will be altered by GSR in a regional and network-specific way.
Furthermore, from the work of (He et al., 2010;Liu et al., 2017;Tong et al., 2018) we know that the GS and vessel BOLD signals (VBS) are highly related, and affect the statistical dependencies in the connectome (Erdo gan et al., 2016). We set out to investigate their unique and joint informative role in the connectome, beyond this high correlation. By applying multiscale partial decomposition of predictive information (Faes et al., 2017) we were able to disambiguate the reduction of uncertainty due to the GS and VBS and found a spatiotemporal modulation across intrinsic connectivity networks and multiple time scales. For the first time, we apply an information theoretical approach to the debate of GS that goes beyond a correlational framework and a nuisance regressor one. This could further help us clarify the effects of GSR, and the information we obtain from the global and vessel BOLD signal at different spatiotemporal scales. Aside from empirical rs-fMRI data, we used other modalities to try and explain the origin of the (predictive) information decomposed pattern of the GS and VBS. In the simulations we tried to explain the spatiotemporal pattern by varying the blood arrival time of simulated rs-fMRI sources. Furthermore, using hemodynamic and calcium mouse recordings we validated the vascular origin of the VBS.
In short, the main goals of this study are: 1) characterizing the effect of GSR on network-level FC estimates in a large dataset; 2) assess the complementary role of global signal and vessels BOLD signal in this modulation; 3) use the framework of (predictive) information decomposition to further examine their joint and unique dynamics. We used a large public dataset of rs-fMRI data, simulations and simultaneous calcium and hemodynamic recordings from mouse data. We would like to emphasize that to the goal is to complement the existing literature and provide an explanation of the GSR-related modulations in the framework of information theory, without explicitly advocating for or against the use of GSR.

HCP rs-fMRI data
We used the minimally cleaned volumetric data from the publicly N. Colenbier et al. NeuroImage 213 (2020) 116699 available 900 Subject release data of the Human Connectome Project (HCP; Van Essen et al., 2012). Subjects underwent neuroimaging protocols, including resting state fMRI, which were collected over two sessions. In this study we considered the data for each subject from the first session, acquired with the left-to-right phase encoding ('rfMRI_R-EST1_LR.nii runs'). All participants gave written consent and experimental procedures were approved by the Institutional Review Board (IRB # 201204036; Title: 'Mapping the Human Connectome: Structure, Function, and Heritability'). Reuse of public data for research purposes is automatically approved by the ethical Committee of the faculty of Psychology and Educational Sciences of Ghent University. We excluded subjects from analysis if the subject exhibited excessive head motion in more than 50% of the frames (details on motion censoring described below). We retained 817 subjects after motion censoring, that were used for further analysis.
Data acquisition and preprocessing. The description of the data acquisition is taken directly from (Barch et al., 2013):"Whole brain EPI acquisitions were acquired with a 32 channel head coil on a modified 3T Siemens skyra with TR ¼ 720 ms, TE ¼ 33.1 ms, flip angle ¼ 52 , BW ¼ 2290 Hz/Px, in-plane FOV ¼ 208 Â 180 mm, 72 slices, 2.0 mm isotropic voxels, with a multi-band acceleration factor of 8. The HCP offers a minimally processed data release. We further performed three preprocessing pipelines to the minimally preprocessed rs-fMRI data to further clean the data: HCP-NO GSR, HCP-GSR and HCP-PID. In the minimally preprocessed pipeline data underwent gradient distortion correction, head motion correction via volume re-alignment, registration to the structural image, bias-field correction and 4D image intensity normalization . We then further (1) de-trended the images; (2) performed censoring; (3) applied nuisance regression (different regressors for the three pipelines); (4) interpolation of censored frames using least squares spectral estimation of the values at censored frames (Power et al., 2014) so that continuous data can be put through a band-pass filter without re-introducing nuisance signals (Hallquist et al., 2013) or contaminating frames near high motion frames (Carp, 2013;Power et al., 2012); (5) band-pass filtering between 0.008 Hz and 0.1 Hz using the fast Fourier Transform rest_bandpass function from the REST toolbox (Song et al., 2011); the HCP-PID pipeline was not filtered see elaboration below); and (6) time series extraction.
Censoring was performed in the following way: framewise displacement (FD; Jenkinson et al., 2012) and root-mean-square of voxel wise differentiated signal (DVARS; Power et al., 2012) were estimated using fsl_motion_outliers. Volumes with FD > 0.20 mm or DVARS >75 were marked as censored frames (outliers). One frame before and two frames after these volumes were also flagged as censored frames. Uncensored segments with fewer than five contiguous volumes were also labeled as censored frames (Gordon et al., 2016). The thresholds for marking censored frames were chosen based on recent studies using the rs-fMRI HCP data (Kong et al., 2018;Li et al., 2019b). We removed subjects from further analysis if the BOLD runs had more than half of the volumes labeled as censored.
Nuisance regression was performed in the following way for each pipeline: in the HCP-NO GSR and HCP-PID pipeline, regressors consisted of: (1) 24 head -motion parameters; (2) white matter signal; and (3) cerebrospinal fluid. In the HCP-GSR pipeline, the GS was calculated and included as an additional regressor at the time of nuisance regression, in order to perform all nuisance regression at the same time. Here we define the GS as the average BOLD signal across all voxels in the brain, excluding the voxels contained in the BOLD vessel mask. In the other pipelines, the HCP-NO GSR and HCP-PID, the GS was calculated during time series extraction. The BOLD vessel mask for each subject is defined as the mask that contains BOLD signal from the vessels. They are created by thresholding the preprocessed T1w/T2w ratio images. The weighted images are used as the contrast of blood vessels is increased in the image by using the ratio of the T1-and T2-weighted images (Glasser and Van Essen, 2011). Fig. 1 depicts the average group mask of all individual BOLD vessel masks. It's clear that on average we extract vessel signal from some major arteries and sinuses. The extracted vessel signal from each subject is on average a mixture from a) arteries such as the Circle of Willis, basilar, cerebral and anterior arteries. B) cerebral arteries supplying the cerebellum with blood. C) signal from sinuses such as the cerebral vein of Galen, the straight sinus and the internal cerebral vein. Important to note is that for all pipelines the censored frames were ignored when computing the regressors (Power et al., 2014).
Importantly, no band-pass filtering was applied on the HCP-PID pipeline after nuisance regression. Here, partial information decomposition (PID) was applied directly after the preprocessing. The PID applied here is a multi-scale extension that computes terms of information transfer at different time scales. It's based on vector autoregressive models and rescales the original time series from the coefficients of the full model. Further filtering is thus not recommended in this case. See the section 'partial information decomposition' below for more details on the method.
Time series were extracted for the three pipelines as follows: The global signal (GS) was calculated as mentioned above. The vessel BOLD signal (VBS) was calculated by averaging the extracted BOLD signal from the vessels, using the BOLD vessels mask. Additionally, signals were averaged in 278 regions of interest (ROIs) using the Shen parcellation (Shen et al., 2013). Each of the ROIs were then assigned to one of the 9 networks (7 cortical networks, plus subcortical and cerebellum regions) as described in Yeo et al., (2011), to further localize the results in ICNs. We chose the Shen parcellation for the following reasons: (1) It is a widely used parcellation of the brain; (2) the nodes are of comparable size, avoiding differential contributions to the GS due to parcel size; (3) It has enough nodes to have complexity at a global and within network level, but also not too many which could cause too much variability across subjects and lose interpretability; (4) it has a clear mapping to the higher order Yeo ICN parcellation.
After preprocessing, connectivity measures were computed from the HCP-NO GSR and HCP-GSR pipelines and partial information decomposition was computed on the HCP-PID pipeline.
Connectivity. FC was calculated using a correlational and mutual information approach. In both approaches the time series of the 278 ROIs obtained from both pipelines in the previous step: (1) HCP-NO GSR and (2) HCP-GSR data were used to calculate FC. The individual connectivity matrices were not thresholded or binarized, and ordered (columns and rows) according to the 9 networks (7 cortical networks, plus subcortical and cerebellum regions) as described in Yeo et al., (2011). In a final step, the connectivity matrices were averaged over subjects.
In the correlational approach, Pearson correlation coefficients (MATLAB command corr) were used to evaluate FC. With the HCP-NO GSR data, FC was calculated between all pairs of ROI time series to obtain a symmetric 278 Â 278 connectivity matrix for every subject. A similar procedure was followed to obtain the 278 Â 278 connectivity matrix for every subject with the HCP-GSR data. Additionally, the Pearson correlation coefficients between the ROI time series and the GS were calculated. To investigate the contribution of VBS to this relationship, partial correlation coefficients (MATLAB command partialcorr) were calculated between the ROI time series and GS conditioned on the VBS. Pearson and partial correlations were calculated only with the uncensored frames.
In the mutual information (MI) approach, MI and conditioned mutual information (CMI) were calculated using a binless rank-based approach based on copulas (Ince et al., 2017). With the HCP-GSR data, MI was calculated between all pairs of ROI time series to obtain a symmetric 278 Â 278 connectivity matrix for every subject. A similar procedure was followed to obtain the 278 Â 278 connectivity matrix for every subject with the HCP-GSR data. Additionally, the MI between the ROI time series and the GS were calculated. To investigate the contribution of VBS to this relationship, CMI was calculated between the ROI time series and GS conditioned on the VBS. The MI measures used here are non-normalized, resulting from the subtraction of two entropy terms. Mutual and conditional mutual information were calculated only with the uncensored frames.

Calcium and hemodynamic mouse recordings
Data from the work of Matsui and colleagues were used, which were provided by the authors upon request (Matsui et al., 2016(Matsui et al., , 2018a. More specifically, a single brain slice from one mouse which contained simultaneous recordings of calcium and hemodynamics was used. This dataset has the advantage of having a high resolution of the brain and vascular structures. Moreover, the data is less affected by movement due to the mouse being tightly head restraint and induced with a light anesthesia. On the other hand, anesthesia has the possible disadvantage of affecting the neurovascular coupling of the neuronal and hemodynamic activity (Matsui et al., 2018b). As a result, the dynamics of the FC might be different compared to recordings from mice who are awake. For a detailed overview of the animal preparation and acquisition of the simultaneous calcium and hemodynamic recordings see (Matsui et al., 2016(Matsui et al., , 2018a.
Data preprocessing. Briefly, the following steps preprocessing were undertaken: all calcium and hemodynamic images were corrected for possible within-scan motion. The images were then further co-registered, spatially down-sampled by a factor of two and high pass filtered (>0.01 Hz). After filtering, the time series were normalized by subtracting the mean and dividing by the standard deviation. For a detailed overview of preprocessing steps see (Matsui et al., 2016(Matsui et al., , 2018a.
We performed additional preprocessing steps resulting in three pipelines for both the hemodynamic and calcium recordings: mouse-PID, mouse-GSR and mouse-NO GSR. The additional steps and differences between the pipelines are: (1) In the mouse-PID pipeline no further preprocessing steps were conducted; (2) For mouse-NO GSR, the data were band-pass filtered; and (3) for mouse-GSR, the global signal was removed using GSR after band-pass filtering. Band-pass filtering was performed between 0.01 and 0.1 Hz. For the mouse-NO GSR and GSR pipelines, the ROI time-series, GS and average vessels time series were calculated after band-pass filtering the data. The GS was calculated by averaging the time series of all pixels of the full brain slice (and is thus not a signal of a volumetric brain as it happens in the human data here considered). In the hemodynamic recordings the VBS was calculated by averaging the time series that were extracted using a vessel mask (Fig. 2). The same masks were used in the calcium recordings to calculate the average vessel time series. Time series were extracted from 38 hand-drawn ROIs as described in (Matsui et al., 2016, 2018a.
After the preprocessing, connectivity measures were computed for the hemodynamic and calcium recordings from the mouse-NO GSR and mouse-GSR pipelines and partial information decomposition was computed on the mouse-PID pipeline.
Connectivity. For the correlation approach, FC was calculated between all pairs of ROI time series to obtain a 38 Â 38 connectivity matrix for the mouse-no GSR data and mouse-GSR data. A similar procedure was followed for the mutual information approach, resulting in a 38 Â 38 connectivity matrix for the no-GSR and GSR data. Pearson correlation coefficients (MATLAB command corr) were used to evaluate FC, while MI was calculated using a binless rank-based approach (Ince et al., 2017).

Simulations
Data generation and preprocessing. The simulations in this work were based on the work from (Erdo gan et al., 2016) in combination with the publicly available simTB toolbox in order to simulate realistic BOLD signals. SimTB provides tools to simulate functionally related rs-fMRI BOLD signals (Erhardt et al., 2012). Here we focus on generating rs-fMRI BOLD sources that have a neuronal contribution and a systemic "sLFO-blood" contribution that arrives with varying temporal delays (blood arrival time) in the rs-fMRI BOLD sources.
In the simulation setting we considered a resting state design and generated rs-fMRI sources over 1200 time points and a TR of 0.72 s (in line with the HCP data). In total time series for 80 rs-fMRI sources were generated for 100 subjects. We started with generating Gaussian time   Matsui et al. (2016Matsui et al. ( , 2018a in green, and the mask for the vessel BOLD signal in red. series of 1200 timepoints for 80 rs-fMRI sources. An additional source was generated that represents the systemic effect of the sLFO's carried with the blood flow. The sLFO blood source was generated by band-pass filtering a random Gaussian noise signal (zero mean, unit standard deviation 1, 1200 timepoints) between 0.008 and 0.1 Hz. A different delayed version of this sLFO-carrying blood source was then added to each of the time series of the 80 rs-fMRI sources with different temporal shifts ranging from 0 to 10 s. This delay was chosen in order to completely cover the physiological time range for blood to complete its cycle from carotid arteries to the jugular veins which takes about 6-9 s (Crandell et al., 1973).
In order to add neuronal contributions to the rs-fMRI sources, we generated 80 neuronal signals and added these to the rs-fMRI sources. The neuronal signals were divided in 4 blocks of 20 functionally related signals. For each related block of signals, spontaneously generated events were generated that occurred with a probability of 0.5 at each TR. These events were then convolved with a canonical hemodynamic response function (HRF; difference of two gamma functions). Finally, these 80 neuronal signals were added with a contribution of 5% to the 80 rs-fMRI sources (same proportion as in Erdo gan et al., 2016). After generating the time series of the rs-fMRI sources partial information decomposition was computed on this data which we refer to as SIM-PID. The GS here is defined as the average signal of the 80 rs-fMRI sources, while the VBS is defined as the signal of the generated sLFO-blood carrying source.

Partial information decomposition
In multivariate systems it makes sense to investigate the joint effect of several driver variables over one (or more) target(s). Using recent developments in information theory, the partial information decomposition (PID) introduced by Williams and Beer (Williams and Beer, 2010), we can quantify these effects in the amount of information being transferred from drivers (GS and VBS) to targets (ROIs) across multiple temporal scales, using the method described in Faes et al. (2017). PID decomposes the joint and unique information being shared between two drivers and a target process. The most relevant measures in this context are the unique information from each of the drivers to the target, and the redundant and synergistic contributions of the two drivers to the target (Fig. 3B).
Differently from the traditional method of information decomposition (interaction information decomposition; IID), the PID accounts for the fact that synergy and redundancy may co-exist as distinctive nonnegative elements of the information transferred by multiple drivers towards a target(s) and allows quantification of these two elements separately. Whereas in the IID a single measure is provided quantifying the net balance between synergy and redundancy.
In a traditional information-theoretic framework, the directed transfer of information between components of interacting processes is assessed by the transfer entropy (TE). The TE from a driver process to a target process quantifies the amount of information contained in the past states of the drivers that can be used to predict the present state of the target beyond the information contained in its own past. In our study, we consider the ROIs of the brain as the target process and the GS and VBS as the driver processes. Note that in the framework of TE "driver" refers to a variable contributing information on the future of another variable ("target") above and beyond what contributed by the past values of the latter. In this case, the drivers (GS,VBS) timeseries reduce the uncertainty on the future values of the target timeseries (ROIs).
The individual TE of information being transferred individually from GS and VBS towards the ROI(s) is then defined as: I (: ; : |: ) is the conditional mutual information, ROI n is the current state of the ROI, and GS À n and VBS À n are the past states of the ROI, GS and VBS. The joint TE, the information that is being transferred from the two drivers to the target can then be defined as: When the sum of the individual TE's is different than the joint TE, the two drivers (GS and VSB) are interacting with each other while transferring information towards the target (ROI). IDD decompose the joint TE and interpret such interactions in terms of synergy and redundancy. An interaction is synergetic if the two drivers transfer more information towards the target when they are considered together than when they are considered individually. While the opposite is true for a redundant interaction. These types of interactions are represented by the sign of the interaction transfer entropy (ITE), which is obtained by decomposing the joint TE using the IDD and is defined as: Where positive values of the ITE denote synergy (the joint TE is greater than the sum of the two individual TEs). While negative values of the ITE denote redundancy (the sum of individual TEs is larger than the joint TE). As mentioned earlier the disadvantage of using the IDD to decompose joint information from two drivers towards a target, is that the ITE is a single measure that quantifies the net balance between the synergy and redundancy making them mutually exclusive (Fig. 3A).
The PID overcomes this deficiency by allowing the joint TE (T GS; VBS→ROI ) to be decomposed as the sum of four separate terms: the unique TE's which are the contributions from the GS or VBS towards the ROI(s) (U GS→ROI and U VBS→ROI ), the synergistic TE (S GS; VBS→ROI ) and the redundant TE (R GS;VBS→ROI ). The PID is then expressed as: As such the PID provides non-negative measures for redundancy and synergy, accounting for the possibility that they may co-exist as separate elements of information modification. Fig. 3B shows the four PID terms of information transfer from the drivers GS and VBS towards the target ROI(s). From Fig. 3B it's evident that the unique TE from each of the drivers towards the target can be obtained, by removing the shared redundant information with the other driver from its own individual TE: The synergistic TE can then be obtained following the PID expression in (5) and substituting equations (6) and (7) into (4): From equations (6)-(8) it follows that the PID terms can be calculated from the joint and individual TE and the redundant information. The joint and individual TE are known information theoretic measures, but a formal definition for the redundant information is needed in order to complete the PID, making each term computable. In this framework the redundancy is defined as the minimum information transferred individually from each driver towards the target. Redundant information is then defined as: The definition of redundant information in equation (9) satisfies the property that the redundant information is independent of the correlation between the driver processes. Following this definition, it's possible to calculate the information transfer from two drivers (GS and VBS) towards a target (ROIs), quantified by the PID in terms of unique, redundant and synergistic information transfer.
Here we perform an extension of the PID calculated over multiple temporal time scales as described in (Faes et al., 2017). In short, this approach considers the multivariate process (GS, VBS and ROIs) as a vector autoregressive (VAR) process, and makes it possible to mathematically derive VAR parameters as a function of the time scale. The multiscale representation of the VAR process for each time scale is obtained by rescaling the original process at time scale τ. Rescaling consists of two subsequent steps, filtering (using a linear finite impulse response low pass filter using a cutoff frequency 1/2 τ ), and then down-sampling (taking one every τ samples) the original process. State space models are then used to represent the VAR parameters after rescaling the original process at time scale τ. In this way, it's possible to compute the information transfer measures of the PID for each τ, by computing the partial variances of the state space models at each τ. For a full description of the implementation of the PID multiscale approach see (Faes et al., 2017).
For the present study, the measures of information decomposition of the PID were computed at a given amount of assigned time scales τ for the multivariate processes (GS, VBS and ROIs). In this case the two drivers are the GS and the VBS, and the targets are the time series of the individual ROIs. The VAR model order, corresponding to the number of samples used to cover the past of the processes, was set according to the Bayesian Information Criterion (BIC) (Schwarz, 1978). The unique, redundant and synergistic information transfer for the human HCP and simulated data were computed for time scales τ ranging from 1 to 21. On the other hand, the first 100 and 200 scales were considered for the mouse hemodynamic and calcium recordings respectively. For each of these modalities, at each timescale τ the information transfer for all oscillations below the given frequency obtained by (10) were considered.
As a result, the information transfer of the PID terms of all the oscillations below: (1) 0.694 Hz at τ ¼ 1, and 0.033 Hz at τ ¼ 21 for the HCP-PID and sim-PID data; (2) 2.5 Hz at τ ¼ 1, and 0.0250 Hz at τ ¼ 100 for the hemodynamic mouse-PID; (3) 5 Hz at τ ¼ 1, and 0.0250 H at τ ¼ 200 for the calcium mouse-PID, were considered in the PID analyses. This means that with the PID at given time scales τ the information transfer for oscillations that fall in between the band-pass filter of most rs-fMRI studies studying brain connectivity are analyzed (~0.008-0.15 Hz). As such interpretations of unique, redundant and synergetic transfer of the multivariate process (GS, VBS and ROIs) that are of interest for the rs-fMRI domain can be made. PID was performed on the data from each modality from the following pipelines: for the human fMRI data, the HCP-PID pipeline data were used (censored frames were ignored for the calculation of the PID terms). While the mouse-PID pipeline was used for the hemodynamic and calcium recordings, and the SIM-PID data for the simulations. . Network-specific effects of GSR on correlational functional connectivity. The top triangle presents the difference histograms between the Pearson_FC without and with GSR for each of the ICN pairs. The top right histogram shows the scale of Pearson_FC differences with a range À0.5 to 1.5 and is the same for each histogram. The black dotted line represents the zero difference point between GSR/NO GSR. The bottom triangle presents the density scatterplots of the Pearson_FC without GSR (x-axis) and with GSR (y-axis) for each of the ICN pairs. The colormap presents the density from low density (purple) to high density (yellow).

Data and code availability
The toolbox for the partial information decomposition method can be found here: https://github.com/danielemarinazzo/multiscale_PID. The toolbox for Gaussian-Copula mutual information can be found here: htt ps://github.com/robince/gcmi. The simTB toolbox we used for simulations can be found here: http://mialab.mrn.org/software/simtb/, and the custom Matlab script for the simulations and analyses reported in this work, as well as the mapping between the Shen parcellation and the YEO ICNs can be found here: https://github.com/compneuro-da/GSR_PID.
The HCP dataset we used can be found at https://www.humanc onnectome.org/. Some of the figures were made using the Gramm package from Pierre Morel, which can be found at https://github.co m/piermorel/gramm and is described in (Morel, 2018). The toolbox REST that contains the function for band-pass filtering can be found at https://www.nitrc.org/projects/rest/.

Results
The Results section is structured in three main parts: (1) rs-fMRI HCP results, where we show the results of (a) the effects of applying GSR on the FC of the connectome; (b) the relationship of the GS with the VBS and the connectome; (c) applying PID, mapping the unique and joint information of the GS and VBS on the connectome. (2) Simulation results, where we explain the observed pattern of the PID results in rs-fMRI data in terms of blood arrival time using simulations. (3) Results of Fig. 6. Network-specific effects of GSR on mutual information functional connectivity. The top triangle presents the difference histograms between the MI_FC without and with GSR for each of the ICN pairs. The top right histogram shows the scale of MI_FC differences with a range À2 to 3 and is the same for each histogram. The black dotted line represents the zero difference point between GSR/NO GSR. The bottom triangle presents the density scatterplots of the MI_FC without GSR (xaxis) and with GSR (y-axis) for each of the ICN pairs. The colormap presents the density from low density (purple) to high density ( hemodynamic and calcium recordings in mouse, where we show (a) the effects of applying GSR on the FC of the connectome; (b) applying PID, mapping the unique and joint information of GS and vessel signals. The hemodynamic and calcium recordings allow us to disambiguate the origin of the PID patterns in terms of physiological (hemodynamic) and neuronal information.

HCP rs-fMRI
GSR. Fig. 4A shows the Pearson_FC among the 278 ROIs averaged across subjects without (bottom triangle) and with GSR (top triangle) in the HCP dataset. Fig. 4B, shows the MI_FC for the same dataset with and without GSR. In Fig. 4C the results of paired t-tests using permutations between the NO GSR and GSR Pearson_FC are shown. It's clear from this figure that after applying GSR the observed reductions and introductions of anticorrelations in Pearson_FC are significantly different from the NO GSR Pearson_FC across the connectome. Furthermore, the introduction of anticorrelations in the connectome happen in a non-uniform way across regions and networks. This is further depicted in Fig. 5, where the scatterplots and the difference histograms between the NO GSR and GSR Pearson_FC are shown for each ICN pair. This figure shows that the shift in Pearson_FC due to applying GSR is different for each ICN pair. As a result, we can derive that GSR has network specific effects across the connectome. The same is true for the MI_FC, where GSR also reduces the FC across ROI pairs as shown by paired t-tests using permutations between the NO GSR and GSR MI_FC in Fig. 4D. As with Pearson_FC, network-specific effects are found, as shown in Fig. 6 by the scatterplots and difference histograms between the NO GSR and GSR MI_FC for each ICN pair.
Two remarks are in order at this point: mutual information is an unsigned quantity, so the difference cannot be expressed in terms of "negative" correlation induced; there is a nonlinear relation between MI and correlation (see also Fig. 3 in Ince et al., 2017), resulting in an enhanced contrast of strong effects with respect to background values. Taken together, GSR reduces the connectivity across ROI pairs, and the intrinsic connectivity networks are affected in non-uniform ways. Our findings align with previous work in the literature (Gotts et al., 2013;Saad et al., 2012). Recent work by Li and colleagues who also studied the effect of GSR in a large sample found similar specific regional and network-related effects introduced by GSR (Li et al., 2019b).   Fig. 7A shows the time series of the GS and VBS from an example subject which are highly related to each other. On average, the GS and VBS are strongly and positively correlated (mean correlation coefficient ¼ 0.87, this is significantly different from 0 (t(816) ¼ 134, p < 0.001). The statistical test was performed after Fisher-z transforming the correlation coefficients. Fig. 7B shows the distribution of the Pearson correlation between the GS and VBS from all subjects. The distribution is skewed to the right, explaining the reported mean 0.87 of the Pearson correlation. The relationship of 88% of the subjects is equal or higher than 0.80, one standard deviation from the mean. We conclude that in general the GS and VBS from HCP subjects are highly related. For interest we also looked into the relationship of the VBS with a respiration derived measure (respiration volume per time), as both influence the FC (see Supplementary Figs. 1 and 2).

Relationship of the global-, vessel BOLD signal and intrinsic connectivity networks
In Fig. 8, the relationship between the GS and the ICN before and after conditioning on the VBS are shown. The GS has an average correlation coefficient of 0.60 and MI of 0.38 with the ICN before controlling for the contribution of VBS. After conditioning on the VBS there is a reduction of the relationship to an average correlation of 0.34 and MI of 0.13. A paired t-test shows that this reduction in correlation (t(816) ¼ 88, p < 0.001) and MI (t(816) ¼ 62, p < 0.001) is significant with large effect sizes for the correlation (Cohen's d ¼ 2.99) and MI (Cohen's d ¼ 2.18).
These results indicate that the GS and VBS are highly related to each other. Additionally, we show that the VBS is an important contributing factor to the relationship between the GS and ICNs, but does not fully eliminate this relationship. Other groups have found similar results showing that the BOLD signal extracted from arteries and veins are highly correlated with the GS (He et al., 2010;Tong et al., 2018). The sLFO's that travel with the blood flow through the brain are present in the GS and the VBS, explaining why they are highly correlated. Despite being highly correlated, our conditioning results on the VBS indicate that both signals have differential and joint effects on the ICNs. However, another framework is needed to further disambiguate the effects of the GS and VBS on the connectome. By applying the PID framework we were able to disambiguate in space and time the differential and joint predictive information for both signals on the different ROIs comprising each ICN.
Partial information decomposition. In Fig. 9 we depict the terms of the PID, applied from the two drivers (GS and VBS) towards the individual ROIs ordered into the ICNs at each of the 21 timescales, interpreted as cutoff frequencies. Kruskal-Wallis tests were conducted to examine network-specific effects of information transfer towards the different ICN for each term of the PID. For each element (unique information two drivers, synergistic and redundant transfer) the time scale for which the maximum transfer value was attained, was used to test for differences between the ICNs. For the unique transfer of the GS (Chi square ¼ 2405.53, p < 0.001, df ¼ 8), unique transfer of the VBS (Chi square ¼ 2361.2, p < 0.001, df ¼ 8), redundant transfer (Chi square ¼ 2377.7, p < 0.001, df ¼ 8) and synergistic transfer (Chi square ¼ 2485.33, p < 0.001, df ¼ 8), significant differences were found among the nine ICNs (VIS, SM, DA, VA, L, FP, DMN, SBC, CB). These tests show that for each of the PID terms, the information transfer towards the ICNs are different at the time scale of maximum information transfer, confirming network-specific effects.
Furthermore, as depicted in Fig. 10, the significant differences (assessed with a permutation based paired t-test) between the unique information from the GS and VBS to the ROIs are different between the ICNs across timescales. The unique GS information is significant at low time scales for all ICNs while for the CB network is also significant at higher time scales (Figs. 9A and 10). The unique VBS information is lower compared to the unique GS information (Fig. 9A and B). The unique VBS information towards the VIS, SM, DA, VA, FP, SBC and CB is significantly higher at higher time scales than the unique GS information (Fig. 9B). However, it is only more significant in the VIS, SM, DA, VA, L, FP, DMN at higher timescales than those of interest in most rs-fMRI studies (Fig. 10). It's clear from the PID that both synergistic and redundant transfer occur at the same time across ICNs and multiple time scales (Fig. 9C and D).
We observe by applying PID that the information terms are modulated across two different domains. (1) Spatial modulation: both the GS and VBS, despite being highly related, affect the ICNs in different ways as evidenced by the difference in the unique information transfers. At the same time, the joint information (dynamics) of both signals is also important as they modulate the ICNs in terms of synergistic and redundant transfer. This spatial modulation is confirmed by the Kruskal-Wallis tests and from the paired t-test between the unique information of both signals (Fig. 10); (2) Temporal modulation: as evidenced by all terms, transfer happens at multiple timescales and patterns of transfer towards ICNs change depending on the timescale. Additionally, the importance of decomposing information at multiple scales is supported by the fact that unique, synergistic, and redundant contributions attain maximum values at time scales >1 (Fig. 9). Here we focused on the spatiotemporal modulation of each PID terms individually. In order to compare PID terms on the same scale, see Supplementary Fig. 3.

Simulations
With the simulations we attempt to explain the spatial and temporal modulations of the information decomposition observed in the empirical rs-fMRI HCP data. We hypothesized that blood arrival time (BAT) might be driving the observed spatiotemporal modulations of the unique, redundant and synergistic information of the GS and VBS across space and time. In the simulations a common 'blood source' carrying physiological information (sLFO's) was mixed with interacting rs-fMRI sources and arrived with varying times at the sources.
Partial Information Decomposition. In Fig. 11, we depict the terms of the PID, applied from the two drivers (GS and VBS) towards the simulated rs-fMRI sources ordered according to BAT at each of the 21 timescales, interpreted as cutoff frequencies. For each of the PID terms depicted in Fig. 11A (unique information GS), Fig. 11B (unique information VBS), Fig. 11C (redundancy) and Fig. 11D (synergy), we observe a spatiotemporal modulation in the information transfer patterns that can be explained by the variation of BAT of the sLFO-blood component in the rs-fMRI sources. Fig. 11A, shows the unique information from the GS towards the sources. We observe that when the blood source arrives early in the rs-fMRI sources information transfer is present in lower amounts across lower and higher timescales, while with increasing BAT the information transfer is observed at lower time scales where it attains its maximum. Sources with late BAT show almost no transfer across low timescales and some transfer at higher timescales. Fig. 11B shows the unique information from the VBS towards the sources. A small amount of information transfer is observed at low timescales with an early BAT. With a linearly increasing BAT, unique VBS information transfer is observed at higher timescales where it attains it maximum. The information transfer decreases again with the BAT reaching its maximum delay at low and high timescales. Note that in the simulations the GS and VBS are present in different amounts across timescales and sources, as is observed in empirical rs-fMRI data ( Fig. 9A  and B).
However, at the same time both signals show joint dynamics of information transfer towards the simulated rs-fMRI sources. This is clear from Fig. 11C and D, which show the redundant and synergistic transfer, respectively. In Fig. 11C, the maximum redundant transfer is attained in sources with an early BAT at low time scales. With an increasing BAT we observe redundant transfer at lower and higher timescales, while the redundant transfer decreases again when the BAT reaches its maximum delay. In Fig. 11D the synergistic transfer of information is shown. With an increasing BAT in the sources, the synergistic transfer is observed at lower and higher timescales while also attaining its maximum. The synergistic transfer then decreases across lower and higher timescales when the BAT reaches its maximum delay. The most left panels of Fig. 11 show the mutual information for the rs-fMRI sources band-pass filtered between 0.008 and 0.1 Hz between the GS and VBS and the sources ordered according to BAT. The MI linearly decreases with the VBS with increasing BAT, while with the GS the MI peaks around 6s BAT and then decreases again. At the same time the information in this plot is complementary to the PID quantities as the unique, redundant, and synergistic information terms quantities (Fig. 11A, B,C,D) are related to the fundamental mutual information measures (Williams and Beer, 2010).
We observe a spatiotemporal modulation of the PID terms. Just as the PID results from the rs-fMRI HCP data, the GS and VBS are present across different timescales and regions (here simulated rs-fMRI sources). With the simulations we tried to explain the spatiotemporal PID pattern observed in Fig. 9, by varying the BAT of a systemic sLFO's source in interacting rs-fMRI sources. The results between the modalities are similar in the sense that we observe a difference in the amount of presence of the unique and joint information of the GS and VBS across networks and time. Here we focused on the spatiotemporal modulation of each PID terms individually. In order to compare PID terms on the same scale, see Supplementary Fig. 4. The patterns are not identical to those of the HCP rs-fMRI data because: (1) here we simulated 80 interacting rs-fMRI sources (similarly to what we would obtain by applying ICA), and not 278 regions that are functionally connected in ICNs; (2) The sources are ordered according to their BAT, while in the rs-fMRI HCP the data is ordered in ICNs and not according to BAT (Fig. 9). Yet, these simulations give a clear indication that the blood arrival time in regions in part explains the observed spatiotemporal pattern in empirical rs-fMRI data.

Simultaneous calcium and hemodynamic recordings
GSR. Fig. 12A depicts the Pearson_FC among the 38 ROIs without (bottom triangle) and with GSR (top triangle) in the single mouse hemodynamic recordings. Fig. 12B shows the MI_FC for the same dataset with and without GSR. Applying GSR reduces Pearson_FC across almost all ROI pairs. Anticorrelations are introduced across ROIs in a nonuniform way across the connectome. The same is true for the MI_FC, where GSR reduces the FC across ROI pairs. These findings are similar to the results we found with the HCP rs-fMRI dataset (Fig. 4). Namely, GSR has local and network-specific effects in reducing FC.
Partial Information Decomposition. In Fig. 13A-D, we depict the results of the PID applied to the hemodynamic recordings. The terms of the PID, applied from the two drivers (hemodynamic GS, hemodynamic VBS) towards the individual hemodynamic ROI signals at each of the 100 timescales, interpreted as cutoff frequencies. For each of the PID components, we observe a spatial and temporal modulation of information transfer towards the individual ROI signals of the mouse. From the unique information from the GS and VBS (Fig. 13A and B) we observe transfer that is different towards the individual ROIs and timescales. Furthermore, both signals show joint transfer of redundant and Fig. 11. Partial information decomposition (PID) in simulated rs-fMRI data. Partial information decomposition (PID) in simulated rs-fMRI data of the information transfer from the GS and VBS towards the individual sources. The PID for the 80 rs-fMRI sources ordered according to blood arrival time, is plotted in function of the time scales 1 to 21 in terms of oscillations analyzed. (A) the unique information from the GS to the individual sources. (B) the unique information from the VBS to the individual sources. (C) the redundant transfer from the GS and VBS to the individual sources. (D) the synergistic transfer from the GS and VBS signal to the individual sources. Left panels: the mutual information between GS and ROIS (blue), VBS and ROIs (red) plotted according to BAT. The signals of the 80 rs-fMRI sources were all band-pass filtered between 0.008 and 0.1 Hz before computing the mutual information.
synergistic information (Fig. 13C and D) which is different for the ROIs across lower and higher timescales.
In Fig. 13 E-H, we depict the results of the PID applied to the calcium recordings. The terms of the PID, applied from the two drivers (calcium GS, calcium vessel signal) towards the individual calcium ROI signals at each of the 200 timescales, interpreted as cutoff frequencies. In this case a spatial and temporal modulation is only observed in the unique information from the GS (Fig. 13E) and redundant information (Fig. 13G). The transfer of unique information from the VBS observed in the hemodynamic recordings (Fig. 13B) disappears in calcium recordings (Fig. 13E), as well as most of its synergy with GS, except from short time scales (Fig. 13H). Calcium recordings are involved in vasodilation changes (Amberg and Navedo, 2013;Behringer, 2017). One could make sense of these results if calcium recordings reflect neuronal activity more prominently in contrast to hemodynamic recordings. Therefore calcium vessel signal might be less prone to containing physiological information in contrast to hemodynamic vessel BOLD signal. If we follow this interpretation then: (1) there is no unique information from calcium vessel signal (Fig. 13E) and synergistic transfer compared ( Fig. 13H) with the hemodynamic recordings (Fig. 13B, D), where the hemodynamic VBS does contain physiological information; (2) The pattern of redundant transfer (Fig. 13G) can then most likely be explained by noise in the calcium vessel signal, which shares some component with the neuronal GS and predicts the information towards the individual ROIs. One possible source of noise could be from the optical imaging method. Some fluctuations could be induced in the signals by the intensity of excitation light which could affect the entire field of view if it exists. The pattern of redundant transfer in the hemodynamic recordings (Fig. 13C) is due to the sharing of physiological information in the GS and VBS.
In Fig. 13I-L, we depict the results of the PID applied to a combination of the hemodynamic and calcium recordings, after down-sampling the latter by a factor 2 to match the former. The terms of the PID, applied from the two drivers (calcium GS, hemodynamic VBS) towards the individual calcium ROI signals at each of the 100 timescales, interpreted as cutoff frequencies. In this case, we observe a spatial and temporal modulation of transfer in the unique information from the GS (Fig. 13I) and redundant and synergistic information ( Fig. 13K and L). Here the effect of unique information from the hemodynamic VBS disappears (Fig. 13J). This shows us that there is no transfer of information of hemodynamic VBS towards calcium ROI signals. This gives us evidence that neuronal functional activity is unaffected by the fluctuations of hemodynamic VBS.
When we consider the PID results using the hemodynamic GS and VBS as drivers, we observe similar findings as with the other modalities: the GS and VBS are present in different amounts across regions and time.
However, the exact PID spatiotemporal modulation pattern of the mouse recordings and the HCP rs-fMRI are different, and this should be expected as: (1) In the mouse data there has been effort to exclude movement using a head restraint and inducing a light anesthesia, while in the HCP data no such efforts have been made to eliminate the effects of movement; (2) We expect differences in the shape of the hemodynamic response function between the species and this influences the decomposed information pattern; (3) There is no clear mapping between the parcellation of the mouse to the Shen parcellation used in HCP rs-fMRI data; (4) The induced anesthesia could have the possible disadvantage of affecting the neurovascular coupling of the neuronal and hemodynamic activity (Matsui et al., 2018b). Results of awake mice might be different and might be more generalizable to human data. Here we focused on the spatiotemporal modulation of each PID terms individually. In order to compare PID terms on the same scale, see Supplementary Fig. 5. Yet, we believe that the results of calcium and hemodynamic mouse recordings from a single slice give us valuable insights. The results confirm that the informational contributions of GS and VBS are present across species, although not in the same pattern for reasons mentioned above, in different amounts across regions and timescales ( Fig. 13A-D). Furthermore, they show that the transfer of information from the VBS is of vascular nature in the hemodynamic recordings, as the information transfer disappears in calcium recordings in which the physiological/ hemodynamic transfer function is not to be expected (Fig. 13F, J).

Discussion
In this paper we learned that, together with their very high correlation, the global signal and the BOLD signals in the vessels provide both unique and shared contributions to the reduction in uncertainty on the dynamics of cortical areas. These differential actions contribute to explaining the effect of GSR on brain-wide intrinsic FC. Additionally, the presence of unique information from both signals could indicate that it is possible, at least theoretically, to separate the two types of signals, hopefully keeping the neural and cognitive correlates, while discarding the physiological nuisance.
Our observations after applying GSR to a large rs-fMRI dataset are consistent with the literature. Using a correlational and mutual information framework, we observed that GSR reduces the functional and mutual information connectivity between ROI pairs ( Fig. 4A and B). This The same 38 ROIs were used as described in (Matsui et al., 2016(Matsui et al., , 2018a. (caption on next page) is in line with previous work (Fox et al., 2009;Murphy et al., 2009;Weissenbacher et al., 2009). Moreover, we observed network-specific effects. The reduction in functional and mutual information connectivity was not uniform across regions and networks (Figs. 5 and 6). This is consistent with previous findings from ( Gotts et al., 2013;Saad et al., 2012), and recent work by (Li et al., 2019b) who also studied the effect of GSR in large human rs-fMRI datasets. These findings were verified in hemodynamic mouse recordings (Fig. 12).
In a next step, we investigated the role of the global signal and vessel BOLD signal in this modulation and further mapped their presence across intrinsic connectivity networks and time using PID. First, we found that the GS and VBS are highly related (Fig. 7), as reported in previous work (Tong et al., 2018). According to Tong and colleagues this finding can be explained by the fact that slow low frequency oscillations (sLFO's) are widespread across the brain travelling with the cerebral blood flow. During the calculation of the GS, neuronal activation which is more local is smoothed/cancelled out leaving the sLFO's as the dominant signal in both the GS and VBS. We further found that the VBS is an important contributor between the relationship of the GS and the intrinsic connectivity networks, but does not fully explain the relationship (Fig. 8). So, despite being highly related they both seem to have some unique relatedness to the intrinsic connectivity networks. By applying PID, we were able to disambiguate the unique and joint dynamics of the GS and VBS by quantifying their presence across networks and time.
We observed in empirical rs-fMRI data, simulations, and hemodynamic mouse recordings that the predictive information from the GS and VBS is present in different amounts across time and space (Figs. 9, 11 and 13). The information of the GS and VBS are present across the connectome in unique ( Fig. 9A and B; 11A, B; 13A, B) and joint ( Fig. 9C and D; 11C, D; 13C, D) ways. Furthermore, we confirmed that: (1) the spatiotemporal modulation of the predictive information in human rs-fMRI ( Fig. 9) can be explained in terms of blood arrival time in the simulations (Fig. 11); and (2) the spatiotemporal modulation is due to the presence of physiological information (sLFO's) in hemodynamic signals in the VBS (Figs. 9B and 13B). As a result, the predictive information of the VBS disappears in calcium recordings due to the lack of physiological information in neuronal signals (Fig. 13F).
Taken together, both signals are present in different amounts and with different timing across regions/networks and one should be aware of this when interpreting FC estimates. We know that the GS is related to physiology (vascular component), and this is due to non-neuronal sLFO's travelling with the blood flow (Tong et al., 2018). For denoising strategies aiming to remove non-neuronal artifacts such as motion or-, respiration, it could be beneficial to control for the sLFO's by correcting for blood arrival time.
The classic GSR approach has been shown to be effective in accounting for global motion and respiratory artifacts compared with other de-noising approaches Lydon-Staley et al., 2019;Parkes et al., 2018;Power et al., 2017a,b;Satterthwaite et al., 2017). However, it does not consider the dynamic passage of global systemic effects (sLFO's) throughout the brain. The dynamic GSR method by (Erdo gan et al., 2016) can effectively deal with the sLFO's by correcting for blood arrival time and reduces the impact of physiological noise. As a result, the FC estimates provide enhanced specificity and reflect more the underlying biology rather than spurious noise. However, it is not known if the dynamic GSR method removes all physiological artifacts. It controls for the vascular effect of the sLFO's, but other physiological effect such as Mayer waves that are related to the sympathetic nervous system oscillations, might be another potential source of confounding global signal oscillations (Julien, 2006).
From the work of Saad and colleagues we know that GSR alters the underlying correlational structure in unpredictable ways depending on the contribution of regions/networks to the GS (Saad et al., 2012). From our observations we have shown that physiology is related to the GS and present in different amounts across networks. Following this view, regressing out the GS might selectively introduce new-mixed physiological artifacts in the FC estimates across networks if not controlled for. As a result, these new mixed physiological artifacts could confound the interpretations made from FC estimates. Being aware of this is important as a wave of recent studies have shown great interest in the relationship between behavior and the architecture of the brain uncovered by FC (Dubois and Adolphs, 2016;Finn et al., 2015;Kong et al., 2018;Li et al., 2019a,b;Rosenberg et al., 2015). Recent work by Li and colleagues shows that GSR strengthens the association between FC and behavioral measures. One could argue that the reliable association between FC could be partly explained by the introduced physiological artifacts after GSR. Future work that focuses on the relationship between behavior and connectivity could investigate if controlling for the physiology (sLFO's) by blood arrival time improves the association between behavior and FC even more.
While global fluctuations have been linked to artifacts such as participant motion, respiration and sLFO's, other studies have linked global fluctuations to neuronal origins. Fluctuations of the rs-fMRI GS have been linked to vigilance (Wong et al., 2013), glucose metabolism (Thompson et al., 2016) and arousal mediated by ascending nuclei Turchi et al., 2018). More recently Gutierrez-Barragan and colleagues have shown that each cycle phase of the GS is the sum of differently overlapping recurring brain states that have a neuronal origin (Gutierrez-Barragan et al., 2019). These brain states that govern the spontaneous brain activity can be captured by CAP-based approaches (Gutierrez-Barragan et al., 2019;Karahano glu & Van De Ville, 2015;Matsui et al., 2016). In addition, work examining global signal topography demonstrates structured information in the GS that can explain individual differences in trait-level cognition and behavior, further suggesting that the signal contains strong cognitive influences (Li et al., 2019a). Aside from the introduction of mixed physiological artifacts, we do not know the additional effects of GSR on recurring brain states with a neuronal origin and how this is reflected in FC estimates.
To conclude, the goal of the current study wasn't to give an overview of the perfect denoising pipeline for rs-fMRI data. We specifically chose to use the minimal preprocessed HCP release opposed to the spatially ICA-FIX cleaned data release, as components classified as noise with this approach can have spatial overlap with signal of interest for us such as the blood vessels (Griffanti et al., 2017(Griffanti et al., , 2014. Since part of the VBS signal (and a not clearly defined one) is removed during ICA-FIX cleaning, it wouldn't make sense to analyze the multivariate process (GS,VBS, ROI) using PID on data with such processing. In addition, the current study is also not meant to take a stance regarding the use of GSR. Rather, the goal is to highlight the finding that the GS and physiology are related and present in different amounts across networks. If one wants to apply de-noising strategies such as GSR, one should be aware of the physiological confound in FC estimates. An alternative de-noising strategy that could be considered is dynamic GSR, which effectively reduces the impact of physiological artifacts in FC patterns (Erdo gan et al., 2016). On the other hand, we know that GSR also impacts task-relevant neuronal information, further complicating the interpretation of FC findings. Alternative de-noising approaches that retain neuronal information and remove global artifacts based on temporal ICA might be promising alternatives to GSR , but see also the comment by Power, 2019).
Limitations. To extract the vessel BOLD signal from the human rs-fMRI data, we created masks that contained vessel information after thresholding the T1w/T2w ratio images. However, we realize that the vessel mask we used for vessel BOLD signal extraction has some limitations. A) the HCP brain mask used for the preprocessed data releases of the HCP explicitly removes the superior sagittal and transverse sinuses from the functional data, removing two of the regions with the strongest VBS signal in the brain, making it harder to extract a clean VBS from the preprocessed data releases of the HCP. However, from the group average vessel mask it's clear we still identify some major vessel structures and extract VBS signal (Fig. 1).b) We realize that due to the low spatial resolution of fMRI our masks can contain voxels that could be classified as grey matter voxels. As a result, the BOLD signal of our vessel BOLD signal mask might be a mixture of grey matter-and vessel BOLD signal effects. Despite the possibility of mixing effects in the vessel BOLD signal that could bias our results, we believe the findings using PID to be robust (Figs. 9 and 13): 1) Using mouse data, we found a spatiotemporal modulation of the predictive information just as in human rs-fMRI data (Figs. 9 and 13). Compared with human rs-fMRI data, the mouse data has the advantage of having high spatial resolution, and reduced movement since the animal is constrained in a rig. Due to the higher resolution we were able to create a mask that only contains vessel information. The extracted VBS in this case is not a mixture of neuronal and vessel information. As we find similar in humans that the GS and VBS are differently present across space and time (including a very high correlation of hemodynamic signals in the vessels and global signal, Fig. 13M, N), we are confident that a possible mixing of the VBS is not confounding the current results. Further, we found that the unique predictive information of the hemodynamic VBS is absent in calcium recordings (Fig. 13E) and was only found with hemodynamic recordings in mouse data (Fig. 13B). This is because calcium vessel signal does not contain physiological (sLFO's) information. 2) In PID, the measures explicitly look for contributions above and beyond what is contained in the target, meaning that the simple coexistence of two phenomena in a time series does not count in the PID. 3) In the human rs-fMRI we used a large sample, averaging the results over hundred of subjects. We could expect that the location of grey matter voxels and neuronal information that might contribute to the VBS mask is inconsistent across subjects. Therefore, the expected neuronal mixing effect might be averaged out in the large group analysis we performed. On the other hand, it is still possible that a neuronal contribution to the VBS is present across subjects and potentially biases the results even in a large sample size.
To summarize, even though what we call VBS in this work (human rs-fMRI) might be a mixture of neuronal and vessel information, we are quite confident that the observed spatiotemporal modulation of the PID in human rs-fMRI is due to the physiological information (sLFO's) that is present in the VBS. Work that aims to validate the results in human rs-fMRI, excluding any neuronal contamination effects, might consider extracting the mask for the VBS from subjects that have angiography data available.

Conclusion
We confirmed that GSR reduces FC estimates between regions and networks in a large empirical fMRI dataset. Furthermore, using PID we show that the GS and sLFO's (physiological artifact) are present in different amounts across different timings, regions and networks. Using simulations, we were able to explain the spatiotemporal modulation of the PID in terms of blood arrival time. Thus, correcting for the sLFO's by taking blood arrival time into account might reduce the introduction of physiological artifacts in FC.